| blob | b71934b6411232d21577d82905f2fe12cc403c2d |
2 !WRF:MODEL_LAYER:DYNAMICS
3 !
5 MODULE module_em
7 USE module_model_constants
9 USE module_advect_em, only: advect_u, advect_v, advect_w, advect_scalar, advect_scalar_pd, advect_scalar_mono, &
10 advect_weno_u, advect_weno_v, advect_weno_w, advect_scalar_weno, advect_scalar_wenopd
12 USE module_big_step_utilities_em, only: grid_config_rec_type, calculate_full, couple_momentum, calc_mu_uv, calc_mu_uv_1, &
13 calc_ww_cp, calc_cq, calc_alt, calc_php, set_tend, rhs_ph, horizontal_pressure_gradient, pg_buoy_w, w_damp, &
14 perturbation_coriolis, coriolis, curvature, horizontal_diffusion, horizontal_diffusion_3dmp, &
15 vertical_diffusion_u, vertical_diffusion_v, vertical_diffusion, vertical_diffusion_3dmp, sixth_order_diffusion, &
16 rk_rayleigh_damp, theta_relaxation, vertical_diffusion_mp, zero_tend, zero_tend2d
18 USE module_ieva_em, only: advect_u_implicit, advect_v_implicit, advect_w_implicit, advect_s_implicit, advect_ph_implicit, &
19 chk_ieva, ww_split, calc_mut_new
20 USE module_state_description, only: param_first_scalar, p_qr, p_qv, p_qc, p_qg, p_qi, p_qs, tiedtkescheme,ntiedtkescheme, heldsuarez, &
21 positivedef, gdscheme, g3scheme, gfscheme, kfetascheme, mskfscheme, monotonic, wenopd_scalar, weno_scalar, weno_mom
23 USE module_damping_em, only: held_suarez_damp
25 USE module_dm
27 USE module_llxy
29 USE module_domain, ONLY : domain, get_ijk_from_grid
31 USE module_configure, ONLY: grid_config_rec_type, model_config_rec, model_to_grid_config_rec
33 CONTAINS
35 !------------------------------------------------------------------------
37 SUBROUTINE rk_step_prep ( config_flags, rk_step, &
38 u, v, w, t, ph, mu, &
39 c1h, c2h, c1f, c2f, moist, &
40 ru, rv, rw, ww, php, alt, &
41 muu, muv, &
42 mub, mut, phb, pb, p, al, alb, &
43 cqu, cqv, cqw, &
44 msfux, msfuy, &
45 msfvx, msfvx_inv, msfvy, &
46 msftx, msfty, &
47 fnm, fnp, dnw, rdx, rdy, &
48 n_moist, &
49 ids, ide, jds, jde, kds, kde, &
50 ims, ime, jms, jme, kms, kme, &
51 its, ite, jts, jte, kts, kte )
53 IMPLICIT NONE
56 ! Input data.
58 TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
60 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
61 ims, ime, jms, jme, kms, kme, &
62 its, ite, jts, jte, kts, kte
64 INTEGER , INTENT(IN ) :: n_moist, rk_step
66 REAL , INTENT(IN ) :: rdx, rdy
68 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
69 INTENT(IN ) :: u, &
70 v, &
71 w, &
72 t, &
73 ph, &
74 phb, &
75 pb, &
76 al, &
77 alb
79 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
80 INTENT( OUT) :: ru, &
81 rv, &
82 rw, &
83 ww, &
84 php, &
85 cqu, &
86 cqv, &
87 cqw, &
88 alt
90 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
91 INTENT(IN ) :: p
96 REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist ), INTENT( IN) :: &
97 moist
99 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msftx, &
100 msfty, &
101 msfux, &
102 msfuy, &
103 msfvx, &
104 msfvx_inv, &
105 msfvy, &
106 mu, &
107 mub
109 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: muu, &
110 muv, &
111 mut
113 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm, fnp, dnw
114 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: c1h, c2h, c1f, c2f
116 integer :: k
119 !<DESCRIPTION>
120 !
121 ! rk_step_prep prepares a number of diagnostic quantities
122 ! in preperation for a Runge-Kutta timestep. subroutines called
123 ! by rk_step_prep calculate
124 !
125 ! (1) total column dry air mass (mut, call to calculate_full)
126 !
127 ! (2) total column dry air mass at u and v points
128 ! (muu, muv, call to calculate_mu_uv)
129 !
130 ! (3) mass-coupled velocities for advection
131 ! (ru, rv, and rw, call to couple_momentum)
132 !
133 ! (4) omega (call to calc_ww_cp)
134 !
135 ! (5) moisture coefficients (cqu, cqv, cqw, call to calc_cq)
136 !
137 ! (6) inverse density (alt, call to calc_alt)
138 !
139 ! (7) geopotential at pressure points (php, call to calc_php)
140 !
141 !</DESCRIPTION>
143 CALL calculate_full( mut, mub, mu, &
144 ids, ide, jds, jde, 1, 2, &
145 ims, ime, jms, jme, 1, 1, &
146 its, ite, jts, jte, 1, 1 )
148 CALL calc_mu_uv ( config_flags, &
149 mu, mub, muu, muv, &
150 ids, ide, jds, jde, kds, kde, &
151 ims, ime, jms, jme, kms, kme, &
152 its, ite, jts, jte, kts, kte )
154 CALL couple_momentum( muu, ru, u, msfuy, &
155 muv, rv, v, msfvx, msfvx_inv, &
156 mut, rw, w, msfty, &
157 c1h, c2h, c1f, c2f, &
158 ids, ide, jds, jde, kds, kde, &
159 ims, ime, jms, jme, kms, kme, &
160 its, ite, jts, jte, kts, kte )
162 ! new call, couples V with mu, also has correct map factors. WCS, 3 june 2001
163 CALL calc_ww_cp ( u, v, mu, mub, c1h, c2h, ww, &
164 rdx, rdy, msftx, msfty, &
165 msfux, msfuy, msfvx, msfvx_inv, &
166 msfvy, dnw, &
167 ids, ide, jds, jde, kds, kde, &
168 ims, ime, jms, jme, kms, kme, &
169 its, ite, jts, jte, kts, kte )
171 CALL calc_cq ( moist, cqu, cqv, cqw, n_moist, &
172 ids, ide, jds, jde, kds, kde, &
173 ims, ime, jms, jme, kms, kme, &
174 its, ite, jts, jte, kts, kte )
176 CALL calc_alt ( alt, al, alb, &
177 ids, ide, jds, jde, kds, kde, &
178 ims, ime, jms, jme, kms, kme, &
179 its, ite, jts, jte, kts, kte )
181 CALL calc_php ( php, ph, phb, &
182 ids, ide, jds, jde, kds, kde, &
183 ims, ime, jms, jme, kms, kme, &
184 its, ite, jts, jte, kts, kte )
186 END SUBROUTINE rk_step_prep
188 !-------------------------------------------------------------------------------
190 SUBROUTINE rk_tendency ( config_flags, rk_step, &
191 ru_tend, rv_tend, rw_tend, ph_tend, t_tend, &
192 ru_tendf, rv_tendf, rw_tendf, ph_tendf, t_tendf, &
193 mu_tend, u_save, v_save, w_save, ph_save, &
194 t_save, mu_save, RTHFTEN, &
195 ru, rv, rw, ww, wwE, wwI, &
196 u, v, w, t, ph, &
197 u_old, v_old, w_old, t_old, ph_old, &
198 h_diabatic, phb,t_init, &
199 mu_old, mu, mut, muu, muv, mub, &
200 c1h, c2h, c1f, c2f, &
201 al, ht, alt, p, pb, php, cqu, cqv, cqw, &
202 u_base, v_base, t_base, qv_base, z_base, &
203 msfux, msfuy, msfvx, msfvx_inv, &
204 msfvy, msftx, msfty, &
205 clat, f, e, sina, cosa, &
206 fnm, fnp, rdn, rdnw, &
207 dt, rdx, rdy, khdif, kvdif, xkmhd, xkhh, &
208 diff_6th_opt, diff_6th_factor, &
209 adv_opt, &
210 dampcoef,zdamp,damp_opt,rad_nudge, &
211 cf1, cf2, cf3, cfn, cfn1, n_moist, &
212 non_hydrostatic, top_lid, &
213 u_frame, v_frame, &
214 ids, ide, jds, jde, kds, kde, &
215 ims, ime, jms, jme, kms, kme, &
216 its, ite, jts, jte, kts, kte, &
217 max_vert_cfl, max_horiz_cfl)
219 IMPLICIT NONE
221 ! Input data.
223 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
225 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
226 ims, ime, jms, jme, kms, kme, &
227 its, ite, jts, jte, kts, kte
229 LOGICAL , INTENT(IN ) :: non_hydrostatic, top_lid
231 INTEGER , INTENT(IN ) :: n_moist, rk_step
233 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
234 INTENT(IN ) :: ru, &
235 rv, &
236 rw, &
237 ww, &
238 u, &
239 v, &
240 w, &
241 t, &
242 ph, &
243 u_old, &
244 v_old, &
245 w_old, &
246 t_old, &
247 ph_old, &
248 phb, &
249 al, &
250 alt, &
251 p, &
252 pb, &
253 php, &
254 cqu, &
255 cqv, &
256 t_init, &
257 xkmhd, &
258 xkhh, &
259 h_diabatic
261 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
262 INTENT(OUT ) :: ru_tend, &
263 rv_tend, &
264 rw_tend, &
265 t_tend, &
266 ph_tend, &
267 RTHFTEN, &
268 u_save, &
269 v_save, &
270 w_save, &
271 ph_save, &
272 t_save
274 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
275 INTENT(INOUT) :: ru_tendf, &
276 rv_tendf, &
277 rw_tendf, &
278 t_tendf, &
279 ph_tendf, &
280 cqw
282 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: mu_tend, &
283 mu_save
285 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
286 msfuy, &
287 msfvx, &
288 msfvx_inv, &
289 msfvy, &
290 msftx, &
291 msfty, &
292 clat, &
293 f, &
294 e, &
295 sina, &
296 cosa, &
297 mu_old, &
298 mu, &
299 mut, &
300 mub, &
301 muu, &
302 muv, &
303 ht
305 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm, &
306 fnp, &
307 rdn, &
308 rdnw, &
309 u_base, &
310 v_base, &
311 t_base, &
312 qv_base, &
313 z_base, &
314 c1h, &
315 c2h, &
316 c1f, &
317 c2f
319 REAL , INTENT(IN ) :: rdx, &
320 rdy, &
321 dt, &
322 u_frame, &
323 v_frame, &
324 khdif, &
325 kvdif
326 INTEGER, INTENT( IN ) :: diff_6th_opt
327 REAL, INTENT( IN ) :: diff_6th_factor
328 INTEGER, INTENT( IN ) :: adv_opt
330 INTEGER, INTENT( IN ) :: damp_opt, rad_nudge
332 REAL, INTENT( IN ) :: zdamp, dampcoef
334 REAL, INTENT( OUT ) :: max_horiz_cfl
335 REAL, INTENT( OUT ) :: max_vert_cfl
337 REAL :: kdift, khdq, kvdq, cfn, cfn1, cf1, cf2, cf3
338 INTEGER :: i,j,k
339 INTEGER :: time_step
340 INTEGER :: rk_order
341 REAL :: dt_step
343 ! IEVA declarations
344 REAL , DIMENSION( ims:ime, kms:kme, jms:jme ) :: wwE, wwI
345 REAL , DIMENSION( ims:ime , jms:jme ) :: mut_old, mut_new
346 REAL , DIMENSION( ims:ime , jms:jme ) :: muu_old, muu_new
347 REAL , DIMENSION( ims:ime , jms:jme ) :: muv_old, muv_new
348 LOGICAL :: ieva
350 !<DESCRIPTION>
351 !
352 ! rk_tendency computes the large-timestep tendency terms in the
353 ! momentum, thermodynamic (theta), and geopotential equations.
354 ! These terms include:
355 !
356 ! (1) advection (for u, v, w, theta - calls to advect_u, advect_v,
357 ! advect_w, and advact_scalar).
358 !
359 ! (2) geopotential equation terms (advection and "gw" - call to rhs_ph).
360 !
361 ! (3) buoyancy term in vertical momentum equation (call to pg_buoy_w).
362 !
363 ! (4) Coriolis and curvature terms in u,v,w momentum equations
364 ! (calls to subroutines coriolis, curvature)
365 !
366 ! (5) 3D diffusion on coordinate surfaces.
367 !
368 !</DESCRIPTION>
370 CALL zero_tend ( ru_tend, &
371 ids, ide, jds, jde, kds, kde, &
372 ims, ime, jms, jme, kms, kme, &
373 its, ite, jts, jte, kts, kte )
375 CALL zero_tend ( rv_tend, &
376 ids, ide, jds, jde, kds, kde, &
377 ims, ime, jms, jme, kms, kme, &
378 its, ite, jts, jte, kts, kte )
380 CALL zero_tend ( rw_tend, &
381 ids, ide, jds, jde, kds, kde, &
382 ims, ime, jms, jme, kms, kme, &
383 its, ite, jts, jte, kts, kte )
385 CALL zero_tend ( t_tend, &
386 ids, ide, jds, jde, kds, kde, &
387 ims, ime, jms, jme, kms, kme, &
388 its, ite, jts, jte, kts, kte )
390 CALL zero_tend ( ph_tend, &
391 ids, ide, jds, jde, kds, kde, &
392 ims, ime, jms, jme, kms, kme, &
393 its, ite, jts, jte, kts, kte )
395 CALL zero_tend ( u_save, &
396 ids, ide, jds, jde, kds, kde, &
397 ims, ime, jms, jme, kms, kme, &
398 its, ite, jts, jte, kts, kte )
400 CALL zero_tend ( v_save, &
401 ids, ide, jds, jde, kds, kde, &
402 ims, ime, jms, jme, kms, kme, &
403 its, ite, jts, jte, kts, kte )
405 CALL zero_tend ( w_save, &
406 ids, ide, jds, jde, kds, kde, &
407 ims, ime, jms, jme, kms, kme, &
408 its, ite, jts, jte, kts, kte )
410 CALL zero_tend ( ph_save, &
411 ids, ide, jds, jde, kds, kde, &
412 ims, ime, jms, jme, kms, kme, &
413 its, ite, jts, jte, kts, kte )
415 CALL zero_tend ( t_save, &
416 ids, ide, jds, jde, kds, kde, &
417 ims, ime, jms, jme, kms, kme, &
418 its, ite, jts, jte, kts, kte )
420 CALL zero_tend2d( mu_tend, &
421 ids, ide, jds, jde, 1, 1, &
422 ims, ime, jms, jme, 1, 1, &
423 its, ite, jts, jte, 1, 1 )
425 CALL zero_tend2d( mu_save, &
426 ids, ide, jds, jde, 1, 1, &
427 ims, ime, jms, jme, 1, 1, &
428 its, ite, jts, jte, 1, 1 )
430 CALL nl_get_time_step ( 1, time_step )
432 rk_order = config_flags%rk_ord
433 dt_step = dt / (rk_order - rk_step + 1) ! needed for calculations using sub-rk step
436 ! Code for splitting vertical velocity into ex/im parts
438 CALL WW_SPLIT(wwE, wwI, &
439 u, v, ww, &
440 mut, rdnw, msfty, &
441 c1f, c2f, &
442 rdx, rdy, msfux, msfuy, &
443 msfvx, msfvy, dt, &
444 config_flags, rk_step, &
445 ids, ide, jds, jde, kds, kde, &
446 ims, ime, jms, jme, kms, kme, &
447 its, ite, jts, jte, kts, kte )
449 ! IEVA ON/OFF
451 ieva = CHK_IEVA( config_flags, rk_step )
453 IF( ieva ) THEN
455 ! Need an estimate of the mut at the original ('n') time level...
457 CALL calculate_full( mut_old, mub, mu_old, &
458 ids, ide, jds, jde, 1, 2, &
459 ims, ime, jms, jme, 1, 1, &
460 its, ite, jts, jte, 1, 1 )
462 ! Estimate new temporary column mass
464 CALL calc_mut_new ( u, v, c1h, c2h, &
465 mut_old, muu, muv, mut_new, &
466 dt, rdx, rdy, msftx, msfty, &
467 msfux, msfuy, msfvx, msfvx_inv, &
468 msfvy, rdnw, &
469 ids, ide, jds, jde, kds, kde, &
470 ims, ime, jms, jme, kms, kme, &
471 its, ite, jts, jte, kts, kte )
473 ! Create staggered mass values for u and v points
475 CALL calc_mu_uv_1 ( config_flags, &
476 mut_old, muu_old, muv_old, &
477 ids, ide, jds, jde, kds, kde, &
478 ims, ime, jms, jme, kms, kme, &
479 its, ite, jts, jte, kts, kte )
481 CALL calc_mu_uv_1 ( config_flags, &
482 mut_new, muu_new, muv_new, &
483 ids, ide, jds, jde, kds, kde, &
484 ims, ime, jms, jme, kms, kme, &
485 its, ite, jts, jte, kts, kte )
487 ENDIF
489 ! Okay do normal advection now....using the wwE array
491 IF( (rk_step == rk_order) .and. ( adv_opt == WENO_MOM ) ) THEN
493 CALL advect_weno_u ( u, u , ru_tend, ru, rv, wwE, &
494 c1h, c2h, &
495 mut, time_step, config_flags, &
496 msfux, msfuy, msfvx, msfvy, &
497 msftx, msfty, &
498 fnm, fnp, rdx, rdy, rdnw, &
499 ids, ide, jds, jde, kds, kde, &
500 ims, ime, jms, jme, kms, kme, &
501 its, ite, jts, jte, kts, kte )
503 ELSE
505 CALL advect_u ( u, u , ru_tend, ru, rv, wwE, &
506 c1h, c2h, &
507 mut, time_step, config_flags, &
508 msfux, msfuy, msfvx, msfvy, &
509 msftx, msfty, &
510 fnm, fnp, rdx, rdy, rdnw, &
511 ids, ide, jds, jde, kds, kde, &
512 ims, ime, jms, jme, kms, kme, &
513 its, ite, jts, jte, kts, kte )
514 ENDIF
516 IF( ieva ) THEN
518 CALL advect_u_implicit ( u, u_old, ru_tend, ru, rv, wwI, &
519 c1h, c2h, &
520 muu_old, muu, muu_new, &
521 config_flags, &
522 msfux, msfuy, msfvx, msfvy, &
523 msftx, msfty, &
524 fnm, fnp, &
525 dt_step, &
526 rdx, rdy, rdnw, &
527 ids, ide, jds, jde, kds, kde, &
528 ims, ime, jms, jme, kms, kme, &
529 its, ite, jts, jte, kts, kte )
530 ENDIF
532 IF( (rk_step == rk_order) .and. ( adv_opt == WENO_MOM ) ) THEN
534 CALL advect_weno_v ( v, v , rv_tend, ru, rv, wwE, &
535 c1h, c2h, &
536 mut, time_step, config_flags, &
537 msfux, msfuy, msfvx, msfvy, &
538 msftx, msfty, &
539 fnm, fnp, rdx, rdy, rdnw, &
540 ids, ide, jds, jde, kds, kde, &
541 ims, ime, jms, jme, kms, kme, &
542 its, ite, jts, jte, kts, kte )
544 ELSE
546 CALL advect_v ( v, v , rv_tend, ru, rv, wwE, &
547 c1h, c2h, &
548 mut, time_step, config_flags, &
549 msfux, msfuy, msfvx, msfvy, &
550 msftx, msfty, &
551 fnm, fnp, rdx, rdy, rdnw, &
552 ids, ide, jds, jde, kds, kde, &
553 ims, ime, jms, jme, kms, kme, &
554 its, ite, jts, jte, kts, kte )
555 ENDIF
557 IF( ieva ) THEN
559 CALL advect_v_implicit ( v, v_old, rv_tend, ru, rv, wwI, &
560 c1h, c2h, &
561 muv_old, muv, muv_new, &
562 config_flags, &
563 msfux, msfuy, msfvx, msfvy, &
564 msftx, msfty, &
565 fnm, fnp, &
566 dt_step, &
567 rdx, rdy, rdnw, &
568 ids, ide, jds, jde, kds, kde, &
569 ims, ime, jms, jme, kms, kme, &
570 its, ite, jts, jte, kts, kte )
571 ENDIF
573 IF (non_hydrostatic) THEN
575 IF( (rk_step == rk_order) .and. ( adv_opt == WENO_MOM ) ) THEN
577 CALL advect_weno_w ( w, w, rw_tend, ru, rv, wwE, &
578 c1h, c2h, &
579 mut, time_step, config_flags, &
580 msfux, msfuy, msfvx, msfvy, &
581 msftx, msfty, &
582 fnm, fnp, rdx, rdy, rdn, &
583 ids, ide, jds, jde, kds, kde, &
584 ims, ime, jms, jme, kms, kme, &
585 its, ite, jts, jte, kts, kte )
587 ELSE
589 CALL advect_w ( w, w, rw_tend, ru, rv, wwE, &
590 c1h, c2h, &
591 mut, time_step, config_flags, &
592 msfux, msfuy, msfvx, msfvy, &
593 msftx, msfty, &
594 fnm, fnp, rdx, rdy, rdn, &
595 ids, ide, jds, jde, kds, kde, &
596 ims, ime, jms, jme, kms, kme, &
597 its, ite, jts, jte, kts, kte )
598 ENDIF
600 ENDIF ! non-hydrostatic
602 ! theta flux divergence
604 ! 11/2016 ERM: Use WENO for theta flux on 3rd RK step if using WENO_SCALAR or WENOPD_SCALAR
605 ! to be consistent with other scalar fluxes
606 IF( ( config_flags%scalar_adv_opt == WENO_SCALAR &
607 .or. config_flags%scalar_adv_opt == WENOPD_SCALAR &
608 .or. config_flags%moist_adv_opt == WENO_SCALAR &
609 .or. config_flags%moist_adv_opt == WENOPD_SCALAR ) &
610 .and. (rk_step == rk_order) ) THEN
612 ! also use weno for monotonic scalar option so that the h_ and z_tendency arrays are not needed
614 CALL advect_scalar_weno ( t, t, t_tend, ru, rv, wwE, &
615 c1h, c2h, mut, time_step, &
616 config_flags, &
617 msfux, msfuy, msfvx, msfvy, &
618 msftx, msfty, fnm, fnp, &
619 rdx, rdy, rdnw, &
620 ids, ide, jds, jde, kds, kde, &
621 ims, ime, jms, jme, kms, kme, &
622 its, ite, jts, jte, kts, kte )
623 ELSE
625 CALL advect_scalar ( t, t, t_tend, ru, rv, wwE, &
626 c1h, c2h, &
627 mut, time_step, config_flags, &
628 msfux, msfuy, msfvx, msfvy, &
629 msftx, msfty, fnm, fnp, &
630 rdx, rdy, rdnw, &
631 ids, ide, jds, jde, kds, kde, &
632 ims, ime, jms, jme, kms, kme, &
633 its, ite, jts, jte, kts, kte )
635 ENDIF
637 IF( ieva ) THEN
639 CALL advect_s_implicit ( t, t_old, t_tend, ru, rv, wwI, &
640 c1h, c2h, &
641 mut_old, mut, mut_new, &
642 config_flags, &
643 msfux, msfuy, msfvx, msfvy, &
644 msftx, msfty, &
645 fnm, fnp, &
646 dt_step, &
647 rdx, rdy, rdnw, &
648 ids, ide, jds, jde, kds, kde, &
649 ims, ime, jms, jme, kms, kme, &
650 its, ite, jts, jte, kts, kte )
652 ENDIF
654 IF ( config_flags%cu_physics == GDSCHEME .OR. &
655 config_flags%cu_physics == GFSCHEME .OR. &
656 config_flags%cu_physics == G3SCHEME .OR. &
657 config_flags%cu_physics == NTIEDTKESCHEME ) THEN ! NTiedtke
659 ! theta advection only:
661 CALL set_tend( RTHFTEN, t_tend, msfty, &
662 ids, ide, jds, jde, kds, kde, &
663 ims, ime, jms, jme, kms, kme, &
664 its, ite, jts, jte, kts, kte )
666 END IF
668 CALL rhs_ph( ph_tend, u, v, wwE, ph, ph, phb, w, &
669 mut, muu, muv, &
670 c1f, c2f, &
671 fnm, fnp, &
672 rdnw, cfn, cfn1, rdx, rdy, &
673 msfux, msfuy, msfvx, &
674 msfvx_inv, msfvy, &
675 msftx, msfty, &
676 non_hydrostatic, &
677 config_flags, &
678 ids, ide, jds, jde, kds, kde, &
679 ims, ime, jms, jme, kms, kme, &
680 its, ite, jts, jte, kts, kte )
682 IF( ieva ) THEN
684 CALL advect_ph_implicit ( ph, ph_old, ph_tend, phb, &
685 ru, rv, wwE, wwI, w, &
686 c1f, c2f, &
687 mut, config_flags, &
688 msfux, msfuy, msfvx, msfvy, &
689 msftx, msfty, &
690 fnm, fnp, &
691 dt_step, &
692 rdx, rdy, rdnw, &
693 ids, ide, jds, jde, kds, kde, &
694 ims, ime, jms, jme, kms, kme, &
695 its, ite, jts, jte, kts, kte )
697 IF (non_hydrostatic) THEN
699 CALL advect_w_implicit ( w, w_old, rw_tend, &
700 ru_tend, rv_tend, ht, wwI, &
701 ph, ph_old, ph_tend, &
702 c1f, c2f, cf1, cf2, cf3, &
703 mut_old, mut, mut_new, &
704 config_flags, &
705 msfux, msfuy, msfvx, msfvy, &
706 msftx, msfty, &
707 fnm, fnp, &
708 dt_step, &
709 rdx, rdy, rdn, &
710 ids, ide, jds, jde, kds, kde, &
711 ims, ime, jms, jme, kms, kme, &
712 its, ite, jts, jte, kts, kte )
713 ENDIF
715 ENDIF ! IEVA
717 CALL horizontal_pressure_gradient( ru_tend,rv_tend, &
718 ph,alt,p,pb,al,php,cqu,cqv, &
719 muu,muv,mu,c1h,c2h,fnm,fnp,rdnw,&
720 cf1,cf2,cf3,cfn,cfn1, &
721 rdx,rdy,msfux,msfuy, &
722 msfvx,msfvy,msftx,msfty, &
723 config_flags, non_hydrostatic, &
724 top_lid, &
725 ids, ide, jds, jde, kds, kde, &
726 ims, ime, jms, jme, kms, kme, &
727 its, ite, jts, jte, kts, kte )
729 IF (non_hydrostatic) THEN
730 CALL pg_buoy_w( rw_tend, p, cqw, mu, mub, &
731 c1f,c2f, &
732 rdnw, rdn, g, msftx, msfty, &
733 ids, ide, jds, jde, kds, kde, &
734 ims, ime, jms, jme, kms, kme, &
735 its, ite, jts, jte, kts, kte )
736 ENDIF
738 CALL w_damp ( rw_tend, max_vert_cfl, &
739 max_horiz_cfl, &
740 u, v, ww, w, mut, c1f, c2f, rdnw, &
741 rdx, rdy, msfux, msfuy, msfvx, &
742 msfvy, dt, config_flags, &
743 ids, ide, jds, jde, kds, kde, &
744 ims, ime, jms, jme, kms, kme, &
745 its, ite, jts, jte, kts, kte )
747 IF(config_flags%pert_coriolis) THEN
749 CALL perturbation_coriolis ( ru, rv, rw, &
750 ru_tend, rv_tend, rw_tend, &
751 config_flags, &
752 u_base, v_base, z_base, &
753 muu, muv, c1h, c2h, phb, ph, &
754 msftx, msfty, msfux, msfuy, &
755 msfvx, msfvy, &
756 f, e, sina, cosa, fnm, fnp, &
757 ids, ide, jds, jde, kds, kde, &
758 ims, ime, jms, jme, kms, kme, &
759 its, ite, jts, jte, kts, kte )
760 ELSE
761 CALL coriolis ( ru, rv, rw, &
762 ru_tend, rv_tend, rw_tend, &
763 config_flags, &
764 msftx, msfty, msfux, msfuy, &
765 msfvx, msfvy, &
766 f, e, sina, cosa, fnm, fnp, &
767 ids, ide, jds, jde, kds, kde, &
768 ims, ime, jms, jme, kms, kme, &
769 its, ite, jts, jte, kts, kte )
771 END IF
773 CALL curvature ( ru, rv, rw, u, v, w, &
774 ru_tend, rv_tend, rw_tend, &
775 config_flags, &
776 msfux, msfuy, msfvx, msfvy, &
777 msftx, msfty, &
778 clat, fnm, fnp, rdx, rdy, &
779 ids, ide, jds, jde, kds, kde, &
780 ims, ime, jms, jme, kms, kme, &
781 its, ite, jts, jte, kts, kte )
783 ! Damping option added for Held-Suarez test (also uses lw option HELDSUAREZ)
785 IF (config_flags%ra_lw_physics == HELDSUAREZ) THEN
786 CALL held_suarez_damp ( ru_tend, rv_tend, &
787 ru,rv,p,pb, &
788 ids, ide, jds, jde, kds, kde, &
789 ims, ime, jms, jme, kms, kme, &
790 its, ite, jts, jte, kts, kte )
791 END IF
793 !**************************************************************
794 !
795 ! Next, the terms that we integrate only with forward-in-time
796 ! (evaluate with time t variables).
797 !
798 !**************************************************************
800 forward_step: IF( rk_step == 1 ) THEN
802 diff_opt1 : IF (config_flags%diff_opt .eq. 1) THEN
804 CALL horizontal_diffusion ('u', u, ru_tendf, mut, &
805 c1h, c2h, config_flags, &
806 msfux, msfuy, msfvx, msfvx_inv, &
807 msfvy,msftx, msfty, &
808 khdif, xkmhd, rdx, rdy, &
809 ids, ide, jds, jde, kds, kde, &
810 ims, ime, jms, jme, kms, kme, &
811 its, ite, jts, jte, kts, kte )
813 CALL horizontal_diffusion ('v', v, rv_tendf, mut, &
814 c1h, c2h, config_flags, &
815 msfux, msfuy, msfvx, msfvx_inv, &
816 msfvy,msftx, msfty, &
817 khdif, xkmhd, rdx, rdy, &
818 ids, ide, jds, jde, kds, kde, &
819 ims, ime, jms, jme, kms, kme, &
820 its, ite, jts, jte, kts, kte )
822 CALL horizontal_diffusion ('w', w, rw_tendf, mut, &
823 c1f, c2f, config_flags, &
824 msfux, msfuy, msfvx, msfvx_inv, &
825 msfvy,msftx, msfty, &
826 khdif, xkmhd, rdx, rdy, &
827 ids, ide, jds, jde, kds, kde, &
828 ims, ime, jms, jme, kms, kme, &
829 its, ite, jts, jte, kts, kte )
831 khdq = 3.*khdif
832 CALL horizontal_diffusion_3dmp ( 'm', t, t_tendf, mut, &
833 c1h, c2h, &
834 config_flags, t_init, &
835 msfux, msfuy, msfvx, msfvx_inv, &
836 msfvy, msftx, msfty, &
837 khdq , xkhh, rdx, rdy, &
838 ids, ide, jds, jde, kds, kde, &
839 ims, ime, jms, jme, kms, kme, &
840 its, ite, jts, jte, kts, kte )
842 pbl_test : IF (config_flags%bl_pbl_physics .eq. 0) THEN
844 CALL vertical_diffusion_u ( u, ru_tendf, config_flags, &
845 u_base, c1h, c2h, &
846 alt, muu, rdn, rdnw, kvdif, &
847 ids, ide, jds, jde, kds, kde, &
848 ims, ime, jms, jme, kms, kme, &
849 its, ite, jts, jte, kts, kte )
851 CALL vertical_diffusion_v ( v, rv_tendf, config_flags, &
852 v_base, c1h, c2h, &
853 alt, muv, rdn, rdnw, kvdif, &
854 ids, ide, jds, jde, kds, kde, &
855 ims, ime, jms, jme, kms, kme, &
856 its, ite, jts, jte, kts, kte )
858 IF (non_hydrostatic) &
859 CALL vertical_diffusion ( 'w', w, rw_tendf, config_flags, &
860 c1f, c2f, &
861 alt, mut, rdn, rdnw, kvdif, &
862 ids, ide, jds, jde, kds, kde, &
863 ims, ime, jms, jme, kms, kme, &
864 its, ite, jts, jte, kts, kte )
866 kvdq = 3.*kvdif
867 CALL vertical_diffusion_3dmp ( t, t_tendf, config_flags, t_init, &
868 c1h, c2h, &
869 alt, mut, rdn, rdnw, kvdq , &
870 ids, ide, jds, jde, kds, kde, &
871 ims, ime, jms, jme, kms, kme, &
872 its, ite, jts, jte, kts, kte )
874 ENDIF pbl_test
876 ! Theta tendency computations.
878 END IF diff_opt1
880 IF ( diff_6th_opt .NE. 0 ) THEN
882 CALL sixth_order_diffusion( 'u', u, ru_tendf, mut, dt, &
883 config_flags, c1h, c2h, &
884 diff_6th_opt, diff_6th_factor, &
885 phb, ph, &
886 rdx, rdy, &
887 msftx, msfty, &
888 msfux, msfuy, &
889 msfvx, msfvy, &
890 ids, ide, jds, jde, kds, kde, &
891 ims, ime, jms, jme, kms, kme, &
892 its, ite, jts, jte, kts, kte )
894 CALL sixth_order_diffusion( 'v', v, rv_tendf, mut, dt, &
895 config_flags, c1h, c2h, &
896 diff_6th_opt, diff_6th_factor, &
897 phb, ph, &
898 rdx, rdy, &
899 msftx, msfty, &
900 msfux, msfuy, &
901 msfvx, msfvy, &
902 ids, ide, jds, jde, kds, kde, &
903 ims, ime, jms, jme, kms, kme, &
904 its, ite, jts, jte, kts, kte )
906 IF (non_hydrostatic) &
907 CALL sixth_order_diffusion( 'w', w, rw_tendf, mut, dt, &
908 config_flags, c1f, c2f, &
909 diff_6th_opt, diff_6th_factor, &
910 phb, ph, &
911 rdx, rdy, &
912 msftx, msfty, &
913 msfux, msfuy, &
914 msfvx, msfvy, &
915 ids, ide, jds, jde, kds, kde, &
916 ims, ime, jms, jme, kms, kme, &
917 its, ite, jts, jte, kts, kte )
919 CALL sixth_order_diffusion( 'm', t, t_tendf, mut, dt, &
920 config_flags, c1h, c2h, &
921 diff_6th_opt, diff_6th_factor, &
922 phb, ph, &
923 rdx, rdy, &
924 msftx, msfty, &
925 msfux, msfuy, &
926 msfvx, msfvy, &
927 ids, ide, jds, jde, kds, kde, &
928 ims, ime, jms, jme, kms, kme, &
929 its, ite, jts, jte, kts, kte )
931 ENDIF
933 IF( damp_opt .eq. 2 ) &
934 CALL rk_rayleigh_damp( ru_tendf, rv_tendf, &
935 rw_tendf, t_tendf, &
936 u, v, w, t, t_init, &
937 c1h, c2h, c1f, c2f, &
938 mut, muu, muv, ph, phb, &
939 u_base, v_base, t_base, z_base, &
940 dampcoef, zdamp, &
941 ids, ide, jds, jde, kds, kde, &
942 ims, ime, jms, jme, kms, kme, &
943 its, ite, jts, jte, kts, kte )
945 IF( rad_nudge .eq. 1 ) &
946 CALL theta_relaxation( t_tendf, t, t_init, &
947 mut, c1h, c2h, ph, phb, &
948 t_base, z_base, &
949 ids, ide, jds, jde, kds, kde, &
950 ims, ime, jms, jme, kms, kme, &
951 its, ite, jts, jte, kts, kte )
953 END IF forward_step
955 END SUBROUTINE rk_tendency
957 !-------------------------------------------------------------------------------
959 SUBROUTINE rk_addtend_dry ( ru_tend, rv_tend, rw_tend, ph_tend, t_tend, &
960 ru_tendf, rv_tendf, rw_tendf, ph_tendf, t_tendf, &
961 u_save, v_save, w_save, ph_save, t_save, &
962 mu_tend, mu_tendf, rk_step, c1, c2, &
963 h_diabatic, mut, msftx, msfty, msfux, msfuy, &
964 msfvx, msfvx_inv, msfvy, &
965 ids,ide, jds,jde, kds,kde, &
966 ims,ime, jms,jme, kms,kme, &
967 ips,ipe, jps,jpe, kps,kpe, &
968 its,ite, jts,jte, kts,kte )
970 IMPLICIT NONE
972 ! Input data.
974 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
975 ims, ime, jms, jme, kms, kme, &
976 ips, ipe, jps, jpe, kps, kpe, &
977 its, ite, jts, jte, kts, kte
978 INTEGER , INTENT(IN ) :: rk_step
980 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(INOUT) :: ru_tend, &
981 rv_tend, &
982 rw_tend, &
983 ph_tend, &
984 t_tend, &
985 ru_tendf, &
986 rv_tendf, &
987 rw_tendf, &
988 ph_tendf, &
989 t_tendf
991 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(INOUT) :: mu_tend, &
992 mu_tendf
994 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(IN ) :: u_save, &
995 v_save, &
996 w_save, &
997 ph_save, &
998 t_save, &
999 h_diabatic
1001 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut, &
1002 msftx, &
1003 msfty, &
1004 msfux, &
1005 msfuy, &
1006 msfvx, &
1007 msfvx_inv, &
1008 msfvy
1010 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: c1, c2
1012 ! Local
1013 INTEGER :: i, j, k
1015 !<DESCRIPTION>
1016 !
1017 ! rk_addtend_dry constructs the full large-timestep tendency terms for
1018 ! momentum (u,v,w), theta and geopotential equations. This is accomplished
1019 ! by combining the physics tendencies (in *tendf; these are computed
1020 ! the first RK substep, held fixed thereafter) with the RK tendencies
1021 ! (in *tend, these include advection, pressure gradient, etc;
1022 ! these change each rk substep). Output is in *tend.
1023 !
1024 !</DESCRIPTION>
1026 ! Finally, add the forward-step tendency to the rk_tendency
1028 ! u/v/w/save contain bc tendency that needs to be multiplied by msf
1029 ! (u by msfuy, v by msfvx)
1030 ! before adding it to physics tendency (*tendf)
1031 ! For momentum we need the final tendency to include an inverse msf
1032 ! physics/bc tendency needs to be divided, advection tendency already has it
1034 ! For scalars we need the final tendency to include an inverse msf (msfty)
1035 ! advection tendency is OK, physics/bc tendency needs to be divided by msf
1037 DO j = jts,MIN(jte,jde-1)
1038 DO k = kts,kte-1
1039 DO i = its,ite
1040 ! multiply by my to uncouple u
1041 IF(rk_step == 1)ru_tendf(i,k,j) = ru_tendf(i,k,j) + u_save(i,k,j)*msfuy(i,j)
1042 ! divide by my to couple u
1043 ru_tend(i,k,j) = ru_tend(i,k,j) + ru_tendf(i,k,j)/msfuy(i,j)
1044 ENDDO
1045 ENDDO
1046 ENDDO
1048 DO j = jts,jte
1049 DO k = kts,kte-1
1050 DO i = its,MIN(ite,ide-1)
1051 ! multiply by mx to uncouple v
1052 IF(rk_step == 1)rv_tendf(i,k,j) = rv_tendf(i,k,j) + v_save(i,k,j)*msfvx(i,j)
1053 ! divide by mx to couple v
1054 rv_tend(i,k,j) = rv_tend(i,k,j) + rv_tendf(i,k,j)*msfvx_inv(i,j)
1055 ENDDO
1056 ENDDO
1057 ENDDO
1059 DO j = jts,MIN(jte,jde-1)
1060 DO k = kts,kte
1061 DO i = its,MIN(ite,ide-1)
1062 ! multiply by my to uncouple w
1063 IF(rk_step == 1)rw_tendf(i,k,j) = rw_tendf(i,k,j) + w_save(i,k,j)*msfty(i,j)
1064 ! divide by my to couple w
1065 rw_tend(i,k,j) = rw_tend(i,k,j) + rw_tendf(i,k,j)/msfty(i,j)
1066 IF(rk_step == 1)ph_tendf(i,k,j) = ph_tendf(i,k,j) + ph_save(i,k,j)
1067 ! divide by my to couple scalar
1068 ph_tend(i,k,j) = ph_tend(i,k,j) + ph_tendf(i,k,j)/msfty(i,j)
1069 ENDDO
1070 ENDDO
1071 ENDDO
1073 DO j = jts,MIN(jte,jde-1)
1074 DO k = kts,kte-1
1075 DO i = its,MIN(ite,ide-1)
1076 IF(rk_step == 1)t_tendf(i,k,j) = t_tendf(i,k,j) + t_save(i,k,j)
1077 ! divide by my to couple theta
1078 t_tend(i,k,j) = t_tend(i,k,j) + t_tendf(i,k,j)/msfty(i,j) &
1079 + (c1(k)*mut(i,j)+c2(k))*h_diabatic(i,k,j)/msfty(i,j)
1080 ! divide by my to couple heating
1081 ENDDO
1082 ENDDO
1083 ENDDO
1085 DO j = jts,MIN(jte,jde-1)
1086 DO i = its,MIN(ite,ide-1)
1087 ! mu tendencies not coupled with 1/msf
1088 MU_TEND(i,j) = MU_TEND(i,j) + MU_TENDF(i,j)
1089 ENDDO
1090 ENDDO
1092 END SUBROUTINE rk_addtend_dry
1094 !-------------------------------------------------------------------------------
1096 SUBROUTINE rk_scalar_tend ( scs, sce, config_flags, &
1097 tenddec, &
1098 rk_step, dt_step, &
1099 ru, rv, ww, wwE, wwI, &
1100 u_old, v_old, &
1101 mut, mub, mu_old, &
1102 c1h, c2h, c1f, c2f, alt, &
1103 scalar_old, scalar, &
1104 scalar_tends, advect_tend, &
1105 h_tendency, z_tendency, &
1106 RQVFTEN, &
1107 base, moist_step, fnm, fnp, &
1108 msfux, msfuy, msfvx, msfvx_inv, &
1109 msfvy, msftx, msfty, &
1110 rdx, rdy, rdn, rdnw, &
1111 khdif, kvdif, xkmhd, &
1112 diff_6th_opt, diff_6th_factor, &
1113 adv_opt, &
1114 phb, ph, &
1115 mix2_off, mix6_off, &
1116 ids, ide, jds, jde, kds, kde, &
1117 ims, ime, jms, jme, kms, kme, &
1118 its, ite, jts, jte, kts, kte )
1120 IMPLICIT NONE
1122 ! Input data.
1124 TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
1126 LOGICAL , INTENT(IN ) :: tenddec ! tendency term
1128 INTEGER , INTENT(IN ) :: rk_step, scs, sce
1129 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1130 ims, ime, jms, jme, kms, kme, &
1131 its, ite, jts, jte, kts, kte
1133 LOGICAL , INTENT(IN ) :: moist_step
1135 REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce ), &
1136 INTENT(IN ) :: scalar, scalar_old
1138 REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce ), &
1139 INTENT(INOUT) :: scalar_tends
1141 REAL, DIMENSION(ims:ime, kms:kme, jms:jme ), INTENT(INOUT) :: advect_tend
1142 REAL, DIMENSION(ims:ime, kms:kme, jms:jme ), INTENT( OUT) :: h_tendency, z_tendency
1144 REAL, DIMENSION(ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RQVFTEN
1146 REAL, DIMENSION(ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: ru, &
1147 rv, &
1148 ww, &
1149 u_old, &
1150 v_old, &
1151 xkmhd, &
1152 alt, &
1153 phb, &
1154 ph
1156 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm, &
1157 fnp, &
1158 rdn, &
1159 rdnw, &
1160 base
1162 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: c1h, c2h
1163 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: c1f, c2f
1165 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
1166 msfuy, &
1167 msfvx, &
1168 msfvx_inv, &
1169 msfvy, &
1170 msftx, &
1171 msfty, &
1172 mub, &
1173 mut, &
1174 mu_old
1176 REAL , INTENT(IN ) :: rdx, &
1177 rdy, &
1178 khdif, &
1179 kvdif
1181 INTEGER, INTENT( IN ) :: diff_6th_opt
1182 REAL, INTENT( IN ) :: diff_6th_factor
1184 REAL, INTENT(IN ) :: dt_step ! This is the local dt for each RK sub-step
1186 INTEGER, INTENT(IN ) :: adv_opt
1187 LOGICAL, INTENT(IN ) :: mix2_off, mix6_off
1189 ! Local data
1191 INTEGER :: im, i,j,k
1192 INTEGER :: time_step
1193 INTEGER :: rk_order
1194 REAL :: dt ! This is the large time step computed below
1196 REAL :: khdq, kvdq, tendency
1198 ! IEVA variables
1199 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: wwE, wwI
1200 REAL, DIMENSION( ims:ime, jms:jme ) :: mut_old
1201 LOGICAL :: ieva
1203 !<DESCRIPTION>
1204 !
1205 ! rk_scalar_tend calls routines that computes scalar tendency from advection
1206 ! and 3D mixing (TKE or fixed eddy viscosities).
1207 !
1208 !</DESCRIPTION>
1210 khdq = khdif/prandtl
1211 kvdq = kvdif/prandtl
1213 rk_order = config_flags%rk_ord
1214 dt = dt_step * (rk_order - rk_step + 1) ! need large time step
1216 ! IEVA
1218 ieva = CHK_IEVA( config_flags, rk_step )
1220 ! Code for splitting vertical velocity into ex/im parts
1222 CALL WW_SPLIT(wwE, wwI, &
1223 u_old, v_old, ww, &
1224 mut, rdnw, msfty, &
1225 c1f, c2f, &
1226 rdx, rdy, msfux, msfuy, &
1227 msfvx, msfvy, dt, &
1228 config_flags, rk_step, &
1229 ids, ide, jds, jde, kds, kde, &
1230 ims, ime, jms, jme, kms, kme, &
1231 its, ite, jts, jte, kts, kte )
1234 IF( ieva ) THEN
1236 ! Need an estimate of the mut at the original ('n') time level...
1238 CALL calculate_full( mut_old, mub, mu_old, &
1239 ids, ide, jds, jde, 1, 2, &
1240 ims, ime, jms, jme, 1, 1, &
1241 its, ite, jts, jte, 1, 1 )
1242 ENDIF
1244 CALL nl_get_time_step ( 1, time_step )
1246 scalar_loop : DO im = scs, sce
1248 CALL zero_tend ( advect_tend(ims,kms,jms), &
1249 ids, ide, jds, jde, kds, kde, &
1250 ims, ime, jms, jme, kms, kme, &
1251 its, ite, jts, jte, kts, kte )
1253 CALL zero_tend ( h_tendency(ims,kms,jms), &
1254 ids, ide, jds, jde, kds, kde, &
1255 ims, ime, jms, jme, kms, kme, &
1256 its, ite, jts, jte, kts, kte )
1258 CALL zero_tend ( z_tendency(ims,kms,jms), &
1259 ids, ide, jds, jde, kds, kde, &
1260 ims, ime, jms, jme, kms, kme, &
1261 its, ite, jts, jte, kts, kte )
1263 CALL nl_get_time_step ( 1, time_step )
1265 IF( (rk_step == rk_order) .and. (adv_opt == POSITIVEDEF) ) THEN
1267 CALL advect_scalar_pd ( scalar(ims,kms,jms,im), &
1268 scalar_old(ims,kms,jms,im), &
1269 advect_tend(ims,kms,jms), &
1270 h_tendency(ims,kms,jms), &
1271 z_tendency(ims,kms,jms), &
1272 ru, rv, wwE, c1h, c2h, &
1273 mut, mub, mu_old, &
1274 time_step, config_flags, tenddec, &
1275 msfux, msfuy, msfvx, msfvy, &
1276 msftx, msfty, fnm, fnp, &
1277 rdx, rdy, rdnw, dt_step, &
1278 ids, ide, jds, jde, kds, kde, &
1279 ims, ime, jms, jme, kms, kme, &
1280 its, ite, jts, jte, kts, kte )
1282 ELSE IF( (rk_step == rk_order) .and. (adv_opt == MONOTONIC) ) THEN
1284 CALL advect_scalar_mono ( scalar(ims,kms,jms,im), &
1285 scalar_old(ims,kms,jms,im), &
1286 advect_tend(ims,kms,jms), &
1287 h_tendency(ims,kms,jms), &
1288 z_tendency(ims,kms,jms), &
1289 ru, rv, wwE, wwI, c1h, c2h, &
1290 mut, mub, mu_old, &
1291 config_flags, tenddec, &
1292 msfux, msfuy, msfvx, msfvy, &
1293 msftx, msfty, fnm, fnp, &
1294 rdx, rdy, rdnw, dt_step, &
1295 ids, ide, jds, jde, kds, kde, &
1296 ims, ime, jms, jme, kms, kme, &
1297 its, ite, jts, jte, kts, kte )
1299 ELSE IF( (rk_step == rk_order) .and. (adv_opt == WENO_SCALAR) ) THEN
1301 CALL advect_scalar_weno ( scalar(ims,kms,jms,im), &
1302 scalar(ims,kms,jms,im), &
1303 advect_tend(ims,kms,jms), &
1304 ru, rv, wwE, c1h, c2h, &
1305 mut, time_step, &
1306 config_flags, &
1307 msfux, msfuy, msfvx, msfvy, &
1308 msftx, msfty, fnm, fnp, &
1309 rdx, rdy, rdnw, &
1310 ids, ide, jds, jde, kds, kde, &
1311 ims, ime, jms, jme, kms, kme, &
1312 its, ite, jts, jte, kts, kte )
1314 ELSEIF( (rk_step == rk_order) .and. (adv_opt == WENOPD_SCALAR) ) THEN
1316 CALL advect_scalar_wenopd ( scalar(ims,kms,jms,im), &
1317 scalar_old(ims,kms,jms,im), &
1318 advect_tend(ims,kms,jms), &
1319 ru, rv, wwE, c1h, c2h, &
1320 mut, mub, mu_old, &
1321 time_step, config_flags, &
1322 msfux, msfuy, msfvx, msfvy, &
1323 msftx, msfty, fnm, fnp, &
1324 rdx, rdy, rdnw, dt_step, & ! dt_step == dt
1325 ids, ide, jds, jde, kds, kde, &
1326 ims, ime, jms, jme, kms, kme, &
1327 its, ite, jts, jte, kts, kte )
1329 ELSE
1331 CALL advect_scalar ( scalar(ims,kms,jms,im), &
1332 scalar(ims,kms,jms,im), &
1333 advect_tend(ims,kms,jms), &
1334 ru, rv, wwE, c1h, c2h, &
1335 mut, time_step, &
1336 config_flags, &
1337 msfux, msfuy, msfvx, msfvy, &
1338 msftx, msfty, fnm, fnp, &
1339 rdx, rdy, rdnw, &
1340 ids, ide, jds, jde, kds, kde, &
1341 ims, ime, jms, jme, kms, kme, &
1342 its, ite, jts, jte, kts, kte )
1344 END IF
1346 IF( ieva ) THEN
1348 CALL advect_s_implicit ( scalar(ims,kms,jms,im), &
1349 scalar_old(ims,kms,jms,im), &
1350 advect_tend(ims,kms,jms), &
1351 ru, rv, wwI, &
1352 c1h, c2h, &
1353 mut_old, mut, mut, &
1354 config_flags, &
1355 msfux, msfuy, msfvx, msfvy, &
1356 msftx, msfty, &
1357 fnm, fnp, &
1358 dt_step, &
1359 rdx, rdy, rdnw, &
1360 ids, ide, jds, jde, kds, kde, &
1361 ims, ime, jms, jme, kms, kme, &
1362 its, ite, jts, jte, kts, kte )
1364 ENDIF
1366 IF((config_flags%cu_physics == GDSCHEME .OR. config_flags%cu_physics == G3SCHEME .OR. &
1367 config_flags%cu_physics == GFSCHEME .OR. &
1368 config_flags%cu_physics == KFETASCHEME .OR. config_flags%cu_physics == MSKFSCHEME .OR. &
1369 config_flags%cu_physics == TIEDTKESCHEME .OR. config_flags%cu_physics == NTIEDTKESCHEME ) & ! Tiedtke
1370 .and. moist_step .and. ( im == P_QV) ) THEN
1372 CALL set_tend( RQVFTEN, advect_tend, msfty, &
1373 ids, ide, jds, jde, kds, kde, &
1374 ims, ime, jms, jme, kms, kme, &
1375 its, ite, jts, jte, kts, kte )
1376 ENDIF
1378 rk_step_1: IF( rk_step == 1 ) THEN
1380 diff_opt1 : IF (config_flags%diff_opt .eq. 1) THEN
1382 IF (.not. mix2_off) &
1383 CALL horizontal_diffusion ( 'm', scalar(ims,kms,jms,im), &
1384 scalar_tends(ims,kms,jms,im), mut, &
1385 c1h, c2h, config_flags, &
1386 msfux, msfuy, msfvx, msfvx_inv, &
1387 msfvy, msftx, msfty, &
1388 khdq , xkmhd, rdx, rdy, &
1389 ids, ide, jds, jde, kds, kde, &
1390 ims, ime, jms, jme, kms, kme, &
1391 its, ite, jts, jte, kts, kte )
1393 pbl_test : IF (config_flags%bl_pbl_physics .eq. 0) THEN
1395 IF( (moist_step) .and. ( im == P_QV)) THEN
1397 CALL vertical_diffusion_mp ( scalar(ims,kms,jms,im), &
1398 scalar_tends(ims,kms,jms,im), &
1399 config_flags, base, c1h, c2h, &
1400 alt, mut, rdn, rdnw, kvdq , &
1401 ids, ide, jds, jde, kds, kde, &
1402 ims, ime, jms, jme, kms, kme, &
1403 its, ite, jts, jte, kts, kte )
1405 ELSE
1407 CALL vertical_diffusion ( 'm', scalar(ims,kms,jms,im), &
1408 scalar_tends(ims,kms,jms,im), &
1409 config_flags, c1h, c2h, &
1410 alt, mut, rdn, rdnw, kvdq, &
1411 ids, ide, jds, jde, kds, kde, &
1412 ims, ime, jms, jme, kms, kme, &
1413 its, ite, jts, jte, kts, kte )
1415 END IF
1417 ENDIF pbl_test
1419 ENDIF diff_opt1
1421 IF ( (diff_6th_opt .NE. 0) .and. (.not. mix6_off) ) THEN
1423 CALL sixth_order_diffusion( 'm', scalar(ims,kms,jms,im), &
1424 scalar_tends(ims,kms,jms,im), &
1425 mut, dt_step, &
1426 config_flags, c1h, c2h, &
1427 diff_6th_opt, diff_6th_factor, &
1428 phb, ph, &
1429 rdx, rdy, &
1430 msftx, msfty, &
1431 msfux, msfuy, &
1432 msfvx, msfvy, &
1433 ids, ide, jds, jde, kds, kde, &
1434 ims, ime, jms, jme, kms, kme, &
1435 its, ite, jts, jte, kts, kte )
1436 ENDIF
1438 ENDIF rk_step_1
1440 END DO scalar_loop
1442 END SUBROUTINE rk_scalar_tend
1444 !-------------------------------------------------------------------------------
1446 SUBROUTINE q_diabatic_add ( scs, sce, &
1447 dt, mut, c1, c2, &
1448 qv_diabatic, qc_diabatic, &
1449 scalar_tends, &
1450 ids, ide, jds, jde, kds, kde, &
1451 ims, ime, jms, jme, kms, kme, &
1452 its, ite, jts, jte, kts, kte )
1454 IMPLICIT NONE
1456 ! Input data.
1458 INTEGER , INTENT(IN ) :: scs, sce
1459 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1460 ims, ime, jms, jme, kms, kme, &
1461 its, ite, jts, jte, kts, kte
1463 REAL, DIMENSION( ims:ime, jms:jme ) , &
1464 INTENT(IN ) :: mut
1466 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: c1, c2
1468 REAL, DIMENSION(ims:ime, kms:kme, jms:jme ), &
1469 INTENT(IN ) :: qv_diabatic, qc_diabatic
1471 REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce ), &
1472 INTENT(INOUT) :: scalar_tends
1474 REAL , INTENT(IN ) :: dt
1476 ! Local data
1478 INTEGER :: im, i,j,k
1480 !<DESCRIPTION>
1481 !
1482 !</DESCRIPTION>
1484 !!! print *,' q_diab add: '
1485 !!! print *,its,MIN(ite,ide-1)
1486 !!! print *,jts,MIN(jte,jde-1)
1487 !!! print *,kts,kte-1
1489 scalar_loop : DO im = scs, sce
1491 IF( im.eq.p_qv )THEN
1492 !!! print *,' qv '
1493 DO j = jts,MIN(jte,jde-1)
1494 DO k = kts,kte-1
1495 DO i = its,MIN(ite,ide-1)
1496 scalar_tends(i,k,j,im) = scalar_tends(i,k,j,im) + qv_diabatic(i,k,j)*(c1(k)*mut(I,J)+c2(k))
1497 ENDDO
1498 ENDDO
1499 ENDDO
1500 ENDIF
1501 IF( im.eq.p_qc )THEN
1502 !!! print *,' qc '
1503 DO j = jts,MIN(jte,jde-1)
1504 DO k = kts,kte-1
1505 DO i = its,MIN(ite,ide-1)
1506 scalar_tends(i,k,j,im) = scalar_tends(i,k,j,im) + qc_diabatic(i,k,j)*(c1(k)*mut(I,J)+c2(k))
1507 ENDDO
1508 ENDDO
1509 ENDDO
1510 ENDIF
1512 END DO scalar_loop
1514 END SUBROUTINE q_diabatic_add
1516 !-------------------------------------------------------------------------------
1518 SUBROUTINE q_diabatic_subtr( scs, sce, &
1519 dt, &
1520 qv_diabatic, qc_diabatic, &
1521 scalar, &
1522 ids, ide, jds, jde, kds, kde, &
1523 ims, ime, jms, jme, kms, kme, &
1524 its, ite, jts, jte, kts, kte )
1526 IMPLICIT NONE
1528 ! Input data.
1530 INTEGER , INTENT(IN ) :: scs, sce
1531 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1532 ims, ime, jms, jme, kms, kme, &
1533 its, ite, jts, jte, kts, kte
1535 REAL, DIMENSION(ims:ime, kms:kme, jms:jme ), &
1536 INTENT(IN ) :: qv_diabatic, qc_diabatic
1538 REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce ), &
1539 INTENT(INOUT) :: scalar
1541 REAL , INTENT(IN ) :: dt
1543 ! Local data
1545 INTEGER :: im, i,j,k
1547 !<DESCRIPTION>
1548 !
1549 !</DESCRIPTION>
1552 !!! print *,' q_diab subtr, dt = : ',dt
1553 !!! print *,its,MIN(ite,ide-1)
1554 !!! print *,jts,MIN(jte,jde-1)
1555 !!! print *,kts,kte-1
1557 scalar_loop : DO im = scs, sce
1559 IF( im.eq.p_qv )THEN
1560 !!! print *,' qv '
1561 DO j = jts,MIN(jte,jde-1)
1562 DO k = kts,kte-1
1563 DO i = its,MIN(ite,ide-1)
1564 scalar(i,k,j,im) = scalar(i,k,j,im) - dt*qv_diabatic(i,k,j)
1565 ENDDO
1566 ENDDO
1567 ENDDO
1568 ENDIF
1569 IF( im.eq.p_qc )THEN
1570 !!! print *,' qc '
1571 DO j = jts,MIN(jte,jde-1)
1572 DO k = kts,kte-1
1573 DO i = its,MIN(ite,ide-1)
1574 scalar(i,k,j,im) = scalar(i,k,j,im) - dt*qc_diabatic(i,k,j)
1575 ENDDO
1576 ENDDO
1577 ENDDO
1578 ENDIF
1580 END DO scalar_loop
1582 END SUBROUTINE q_diabatic_subtr
1585 !-------------------------------------------------------------------------------
1587 SUBROUTINE rk_update_scalar( scs, sce, &
1588 scalar_1, scalar_2, sc_tend, &
1589 advh_t, advz_t, &
1590 advect_tend, &
1591 h_tendency, z_tendency, &
1592 msftx, msfty, c1, c2, &
1593 mu_old, mu_new, mu_base, &
1594 rk_step, dt, spec_zone, &
1595 config_flags, &
1596 tenddec, &
1597 ids, ide, jds, jde, kds, kde, &
1598 ims, ime, jms, jme, kms, kme, &
1599 its, ite, jts, jte, kts, kte )
1601 IMPLICIT NONE
1603 ! Input data.
1605 TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
1607 LOGICAL , INTENT(IN ) :: tenddec
1609 INTEGER , INTENT(IN ) :: scs, sce, rk_step, spec_zone
1610 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1611 ims, ime, jms, jme, kms, kme, &
1612 its, ite, jts, jte, kts, kte
1614 REAL, INTENT(IN ) :: dt
1616 REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce), &
1617 INTENT(INOUT) :: scalar_1, &
1618 scalar_2
1620 REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce), &
1621 INTENT(IN) :: sc_tend
1623 REAL, DIMENSION(ims:ime, kms:kme, jms:jme ), &
1624 INTENT(IN) :: advect_tend
1626 REAL, DIMENSION(ims:ime, kms:kme, jms:jme ), INTENT(INOUT) , OPTIONAL :: advh_t, advz_t ! accumulating for output
1627 REAL, DIMENSION(ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: h_tendency, z_tendency ! from rk_scalar_tend
1629 REAL, DIMENSION(ims:ime, jms:jme ), INTENT(IN ) :: mu_old, &
1630 mu_new, &
1631 mu_base, &
1632 msftx, &
1633 msfty
1635 REAL, DIMENSION(kms:kme ), INTENT(IN ) :: c1, c2
1637 INTEGER :: i,j,k,im
1638 REAL :: sc_middle, msfsq
1639 REAL, DIMENSION(its:ite) :: muold, munew
1641 REAL, DIMENSION(its:ite, kts:kte, jts:jte ) :: tendency
1643 INTEGER :: i_start,i_end,j_start,j_end,k_start,k_end
1644 INTEGER :: i_start_spc,i_end_spc,j_start_spc,j_end_spc,k_start_spc,k_end_spc
1646 !<DESCRIPTION>
1647 !
1648 ! rk_scalar_update advances the scalar equation given the time t value
1649 ! of the scalar and the scalar tendency.
1650 !
1651 !</DESCRIPTION>
1654 !
1655 ! set loop limits.
1657 i_start = its
1658 i_end = min(ite,ide-1)
1659 j_start = jts
1660 j_end = min(jte,jde-1)
1661 k_start = kts
1662 k_end = kte-1
1664 i_start_spc = i_start
1665 i_end_spc = i_end
1666 j_start_spc = j_start
1667 j_end_spc = j_end
1668 k_start_spc = k_start
1669 k_end_spc = k_end
1671 IF( config_flags%nested .or. config_flags%specified ) THEN
1672 IF( .NOT. config_flags%periodic_x)i_start = max( its,ids+spec_zone )
1673 IF( .NOT. config_flags%periodic_x)i_end = min( ite,ide-spec_zone-1 )
1674 j_start = max( jts,jds+spec_zone )
1675 j_end = min( jte,jde-spec_zone-1 )
1676 k_start = kts
1677 k_end = min( kte, kde-1 )
1678 ENDIF
1680 IF ( rk_step == 1 ) THEN
1682 ! replace t-dt values (in scalar_1) with t values scalar_2,
1683 ! then compute new values by adding tendency to values at t
1685 DO im = scs,sce
1687 DO j = jts, min(jte,jde-1)
1688 DO k = kts, min(kte,kde-1)
1689 DO i = its, min(ite,ide-1)
1690 tendency(i,k,j) = 0.
1691 ENDDO
1692 ENDDO
1693 ENDDO
1695 DO j = j_start,j_end
1696 DO k = k_start,k_end
1697 DO i = i_start,i_end
1698 ! scalar was coupled with my
1699 tendency(i,k,j) = advect_tend(i,k,j) * msfty(i,j)
1700 ENDDO
1701 ENDDO
1702 ENDDO
1704 DO j = j_start_spc,j_end_spc
1705 DO k = k_start_spc,k_end_spc
1706 DO i = i_start_spc,i_end_spc
1707 tendency(i,k,j) = tendency(i,k,j) + sc_tend(i,k,j,im)
1708 ENDDO
1709 ENDDO
1710 ENDDO
1712 DO j = jts, min(jte,jde-1)
1714 DO i = its, min(ite,ide-1)
1715 MUOLD(i) = MU_OLD(i,j) + MU_BASE(i,j)
1716 MUNEW(i) = MU_NEW(i,j) + MU_BASE(i,j)
1717 ENDDO
1719 DO k = kts, min(kte,kde-1)
1720 DO i = its, min(ite,ide-1)
1722 scalar_1(i,k,j,im) = scalar_2(i,k,j,im)
1723 scalar_2(i,k,j,im) = ((c1(k)*muold(i)+c2(k))*scalar_1(i,k,j,im) &
1724 + dt*tendency(i,k,j))/(c1(k)*munew(i)+c2(k))
1726 ENDDO
1727 ENDDO
1728 ENDDO
1730 ENDDO
1732 ELSE
1734 ! just compute new values, scalar_1 already at time t.
1736 DO im = scs, sce
1738 DO j = jts, min(jte,jde-1)
1739 DO k = kts, min(kte,kde-1)
1740 DO i = its, min(ite,ide-1)
1741 tendency(i,k,j) = 0.
1742 ENDDO
1743 ENDDO
1744 ENDDO
1746 DO j = j_start,j_end
1747 DO k = k_start,k_end
1748 DO i = i_start,i_end
1749 ! scalar was coupled with my
1750 tendency(i,k,j) = advect_tend(i,k,j) * msfty(i,j)
1751 ENDDO
1752 ENDDO
1753 ENDDO
1755 DO j = j_start_spc,j_end_spc
1756 DO k = k_start_spc,k_end_spc
1757 DO i = i_start_spc,i_end_spc
1758 tendency(i,k,j) = tendency(i,k,j) + sc_tend(i,k,j,im)
1759 ENDDO
1760 ENDDO
1761 ENDDO
1763 DO j = jts, min(jte,jde-1)
1765 DO i = its, min(ite,ide-1)
1766 MUOLD(i) = MU_OLD(i,j) + MU_BASE(i,j)
1767 MUNEW(i) = MU_NEW(i,j) + MU_BASE(i,j)
1768 ENDDO
1770 DO k = kts, min(kte,kde-1)
1771 DO i = its, min(ite,ide-1)
1773 scalar_2(i,k,j,im) = ((c1(k)*muold(i)+c2(k))*scalar_1(i,k,j,im) &
1774 + dt*tendency(i,k,j))/(c1(k)*munew(i)+c2(k))
1776 ENDDO
1777 ENDDO
1779 ! This is separated from the k/i-loop above for better performance
1780 IF ( PRESENT(advh_t) .AND. PRESENT(advz_t) ) THEN
1781 IF(tenddec.and.rk_step.eq.config_flags%rk_ord) THEN
1782 DO k = kts, min(kte,kde-1)
1783 DO i = its, min(ite,ide-1)
1785 advh_t(i,k,j) = advh_t(i,k,j) + (dt*h_tendency(i,k,j)* msfty(i,j))/(c1(k)*munew(i)+c2(k))
1786 advz_t(i,k,j) = advz_t(i,k,j) + (dt*z_tendency(i,k,j)* msfty(i,j))/(c1(k)*munew(i)+c2(k))
1788 ENDDO
1789 ENDDO
1790 END IF
1791 END IF
1793 ENDDO
1795 ENDDO
1797 END IF
1799 END SUBROUTINE rk_update_scalar
1801 !-------------------------------------------------------------------------------
1803 SUBROUTINE rk_update_scalar_pd( scs, sce, &
1804 scalar, sc_tend, &
1805 c1, c2, &
1806 mu_old, mu_new, mu_base, &
1807 rk_step, dt, spec_zone, &
1808 config_flags, &
1809 ids, ide, jds, jde, kds, kde, &
1810 ims, ime, jms, jme, kms, kme, &
1811 its, ite, jts, jte, kts, kte )
1813 IMPLICIT NONE
1815 ! Input data.
1817 TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
1819 INTEGER , INTENT(IN ) :: scs, sce, rk_step, spec_zone
1820 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1821 ims, ime, jms, jme, kms, kme, &
1822 its, ite, jts, jte, kts, kte
1824 REAL, INTENT(IN ) :: dt
1826 REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce), &
1827 INTENT(INOUT) :: scalar, &
1828 sc_tend
1830 REAL, DIMENSION(ims:ime, jms:jme ), INTENT(IN ) :: mu_old, &
1831 mu_new, &
1832 mu_base
1834 REAL, DIMENSION(kms:kme ), INTENT(IN ) :: c1, c2
1836 INTEGER :: i,j,k,im
1837 REAL :: sc_middle, msfsq
1838 REAL, DIMENSION(its:ite) :: muold, munew
1840 REAL, DIMENSION(its:ite, kts:kte, jts:jte ) :: tendency
1842 INTEGER :: i_start,i_end,j_start,j_end,k_start,k_end
1843 INTEGER :: i_start_spc,i_end_spc,j_start_spc,j_end_spc,k_start_spc,k_end_spc
1845 !<DESCRIPTION>
1846 !
1847 ! rk_scalar_update advances the scalar equation given the time t value
1848 ! of the scalar and the scalar tendency.
1849 !
1850 !</DESCRIPTION>
1853 !
1854 ! set loop limits.
1856 i_start = its
1857 i_end = min(ite,ide-1)
1858 j_start = jts
1859 j_end = min(jte,jde-1)
1860 k_start = kts
1861 k_end = kte-1
1863 i_start_spc = i_start
1864 i_end_spc = i_end
1865 j_start_spc = j_start
1866 j_end_spc = j_end
1867 k_start_spc = k_start
1868 k_end_spc = k_end
1870 IF( config_flags%nested .or. config_flags%specified ) THEN
1871 IF( .NOT. config_flags%periodic_x)i_start = max( its,ids+spec_zone )
1872 IF( .NOT. config_flags%periodic_x)i_end = min( ite,ide-spec_zone-1 )
1873 j_start = max( jts,jds+spec_zone )
1874 j_end = min( jte,jde-spec_zone-1 )
1875 k_start = kts
1876 k_end = min( kte, kde-1 )
1877 ENDIF
1879 DO im = scs, sce
1881 DO j = jts, min(jte,jde-1)
1882 DO k = kts, min(kte,kde-1)
1883 DO i = its, min(ite,ide-1)
1884 tendency(i,k,j) = 0.
1885 ENDDO
1886 ENDDO
1887 ENDDO
1889 DO j = j_start_spc,j_end_spc
1890 DO k = k_start_spc,k_end_spc
1891 DO i = i_start_spc,i_end_spc
1892 tendency(i,k,j) = tendency(i,k,j) + sc_tend(i,k,j,im)
1893 sc_tend(i,k,j,im) = 0.
1894 ENDDO
1895 ENDDO
1896 ENDDO
1898 DO j = jts, min(jte,jde-1)
1900 DO i = its, min(ite,ide-1)
1901 MUOLD(i) = MU_OLD(i,j) + MU_BASE(i,j)
1902 MUNEW(i) = MU_NEW(i,j) + MU_BASE(i,j)
1903 ENDDO
1905 DO k = kts, min(kte,kde-1)
1906 DO i = its, min(ite,ide-1)
1908 scalar(i,k,j,im) = ((c1(k)*muold(i)+c2(k))*scalar(i,k,j,im) &
1909 + dt*tendency(i,k,j))/(c1(k)*munew(i)+c2(k))
1910 ENDDO
1911 ENDDO
1912 ENDDO
1914 ENDDO
1916 END SUBROUTINE rk_update_scalar_pd
1918 !------------------------------------------------------------
1920 SUBROUTINE init_zero_tendency(ru_tendf, rv_tendf, rw_tendf, ph_tendf, &
1921 t_tendf, tke_tendf, mu_tendf, &
1922 moist_tendf,chem_tendf,scalar_tendf, &
1923 tracer_tendf,n_tracer, &
1924 n_moist,n_chem,n_scalar,rk_step, &
1925 ids, ide, jds, jde, kds, kde, &
1926 ims, ime, jms, jme, kms, kme, &
1927 its, ite, jts, jte, kts, kte )
1928 !-----------------------------------------------------------------------
1929 IMPLICIT NONE
1930 !-----------------------------------------------------------------------
1932 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1933 ims, ime, jms, jme, kms, kme, &
1934 its, ite, jts, jte, kts, kte
1936 INTEGER , INTENT(IN ) :: n_moist,n_chem,n_scalar,n_tracer,rk_step
1938 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(INOUT) :: &
1939 ru_tendf, &
1940 rv_tendf, &
1941 rw_tendf, &
1942 ph_tendf, &
1943 t_tendf, &
1944 tke_tendf
1946 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(INOUT) :: mu_tendf
1948 REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_moist),INTENT(INOUT)::&
1949 moist_tendf
1951 REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_chem ),INTENT(INOUT)::&
1952 chem_tendf
1953 REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_tracer ),INTENT(INOUT)::&
1954 tracer_tendf
1956 REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_scalar ),INTENT(INOUT)::&
1957 scalar_tendf
1959 ! LOCAL VARS
1961 INTEGER :: im, ic, is
1963 !<DESCRIPTION>
1964 !
1965 ! init_zero_tendency
1966 ! sets tendency arrays to zero for all prognostic variables.
1967 !
1968 !</DESCRIPTION>
1971 CALL zero_tend ( ru_tendf, &
1972 ids, ide, jds, jde, kds, kde, &
1973 ims, ime, jms, jme, kms, kme, &
1974 its, ite, jts, jte, kts, kte )
1976 CALL zero_tend ( rv_tendf, &
1977 ids, ide, jds, jde, kds, kde, &
1978 ims, ime, jms, jme, kms, kme, &
1979 its, ite, jts, jte, kts, kte )
1981 CALL zero_tend ( rw_tendf, &
1982 ids, ide, jds, jde, kds, kde, &
1983 ims, ime, jms, jme, kms, kme, &
1984 its, ite, jts, jte, kts, kte )
1986 CALL zero_tend ( ph_tendf, &
1987 ids, ide, jds, jde, kds, kde, &
1988 ims, ime, jms, jme, kms, kme, &
1989 its, ite, jts, jte, kts, kte )
1991 CALL zero_tend ( t_tendf, &
1992 ids, ide, jds, jde, kds, kde, &
1993 ims, ime, jms, jme, kms, kme, &
1994 its, ite, jts, jte, kts, kte )
1996 CALL zero_tend ( tke_tendf, &
1997 ids, ide, jds, jde, kds, kde, &
1998 ims, ime, jms, jme, kms, kme, &
1999 its, ite, jts, jte, kts, kte )
2001 CALL zero_tend2d( mu_tendf, &
2002 ids, ide, jds, jde, kds, kds, &
2003 ims, ime, jms, jme, kms, kms, &
2004 its, ite, jts, jte, kts, kts )
2006 ! DO im=PARAM_FIRST_SCALAR,n_moist
2007 DO im=1,n_moist ! make sure first one is zero too
2008 CALL zero_tend ( moist_tendf(ims,kms,jms,im), &
2009 ids, ide, jds, jde, kds, kde, &
2010 ims, ime, jms, jme, kms, kme, &
2011 its, ite, jts, jte, kts, kte )
2012 ENDDO
2014 ! DO ic=PARAM_FIRST_SCALAR,n_chem
2015 DO ic=1,n_chem ! make sure first one is zero too
2016 CALL zero_tend ( chem_tendf(ims,kms,jms,ic), &
2017 ids, ide, jds, jde, kds, kde, &
2018 ims, ime, jms, jme, kms, kme, &
2019 its, ite, jts, jte, kts, kte )
2020 ENDDO
2022 ! DO ic=PARAM_FIRST_SCALAR,n_tracer
2023 DO ic=1,n_tracer ! make sure first one is zero too
2024 CALL zero_tend ( tracer_tendf(ims,kms,jms,ic), &
2025 ids, ide, jds, jde, kds, kde, &
2026 ims, ime, jms, jme, kms, kme, &
2027 its, ite, jts, jte, kts, kte )
2028 ENDDO
2030 ! DO ic=PARAM_FIRST_SCALAR,n_scalar
2031 DO ic=1,n_scalar ! make sure first one is zero too
2032 CALL zero_tend ( scalar_tendf(ims,kms,jms,ic), &
2033 ids, ide, jds, jde, kds, kde, &
2034 ims, ime, jms, jme, kms, kme, &
2035 its, ite, jts, jte, kts, kte )
2036 ENDDO
2038 END SUBROUTINE init_zero_tendency
2040 !===================================================================
2043 SUBROUTINE dump_data( a, field, io_unit, &
2044 ims, ime, jms, jme, kms, kme, &
2045 ids, ide, jds, jde, kds, kde )
2046 implicit none
2047 integer :: ims, ime, jms, jme, kms, kme, &
2048 ids, ide, jds, jde, kds, kde
2049 real, dimension(ims:ime, kms:kme, jds:jde) :: a
2050 character :: field
2051 integer :: io_unit
2053 integer :: is,ie,js,je,ks,ke
2055 !<DESCRIPTION
2056 !
2057 ! quick and dirty debug io utility
2058 !
2059 !</DESCRIPTION
2061 is = ids
2062 ie = ide-1
2063 js = jds
2064 je = jde-1
2065 ks = kds
2066 ke = kde-1
2068 if(field == 'u') ie = ide
2069 if(field == 'v') je = jde
2070 if(field == 'w') ke = kde
2072 write(io_unit) is,ie,ks,ke,js,je
2073 write(io_unit) a(is:ie, ks:ke, js:je)
2075 end subroutine dump_data
2077 !-----------------------------------------------------------------------
2079 SUBROUTINE calculate_phy_tend (config_flags,c1,c2, &
2080 mut,muu,muv,pi3d, &
2081 RTHRATEN, &
2082 RUBLTEN,RVBLTEN,RTHBLTEN, &
2083 RQVBLTEN,RQCBLTEN,RQIBLTEN, &
2084 RUCUTEN,RVCUTEN,RTHCUTEN, &
2085 RQVCUTEN,RQCCUTEN,RQRCUTEN, &
2086 RQICUTEN,RQSCUTEN, &
2087 RUSHTEN,RVSHTEN,RTHSHTEN, &
2088 RQVSHTEN,RQCSHTEN,RQRSHTEN, &
2089 RQISHTEN,RQSSHTEN,RQGSHTEN, &
2090 RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RQVNDGDTEN, &
2091 RMUNDGDTEN, &
2092 scalar, scalar_tend, num_scalar, &
2093 tracer, tracer_tend, num_tracer, &
2094 ids,ide, jds,jde, kds,kde, &
2095 ims,ime, jms,jme, kms,kme, &
2096 its,ite, jts,jte, kts,kte )
2097 !-----------------------------------------------------------------------
2098 IMPLICIT NONE
2100 TYPE(grid_config_rec_type), INTENT(IN) :: config_flags
2102 INTEGER, INTENT(IN ) :: num_scalar, num_tracer
2103 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
2104 ims,ime, jms,jme, kms,kme, &
2105 its,ite, jts,jte, kts,kte
2107 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
2108 INTENT(IN ) :: pi3d
2110 REAL, DIMENSION( kms:kme ) , &
2111 INTENT(IN ) :: c1, c2
2113 REAL, DIMENSION( ims:ime, jms:jme ) , &
2114 INTENT(IN ) :: mut, &
2115 muu, &
2116 muv
2119 ! radiation
2121 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
2122 INTENT(INOUT) :: RTHRATEN
2124 ! cumulus
2126 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
2127 INTENT(INOUT) :: &
2128 RUCUTEN, &
2129 RVCUTEN, &
2130 RTHCUTEN, &
2131 RQVCUTEN, &
2132 RQCCUTEN, &
2133 RQRCUTEN, &
2134 RQICUTEN, &
2135 RQSCUTEN, &
2136 RUSHTEN, &
2137 RVSHTEN, &
2138 RTHSHTEN, &
2139 RQVSHTEN, &
2140 RQCSHTEN, &
2141 RQRSHTEN, &
2142 RQISHTEN, &
2143 RQSSHTEN, &
2144 RQGSHTEN
2146 ! pbl
2148 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
2149 INTENT(INOUT) :: RUBLTEN, &
2150 RVBLTEN, &
2151 RTHBLTEN, &
2152 RQVBLTEN, &
2153 RQCBLTEN, &
2154 RQIBLTEN
2156 ! fdda
2158 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
2159 INTENT(INOUT) :: RUNDGDTEN, &
2160 RVNDGDTEN, &
2161 RTHNDGDTEN, &
2162 RQVNDGDTEN
2163 REAL, DIMENSION( ims:ime, jms:jme ) , &
2164 INTENT(INOUT) :: RMUNDGDTEN
2166 ! 4d arrays
2168 REAL ,DIMENSION(ims:ime,kms:kme,jms:jme,num_tracer),INTENT(INOUT) :: tracer
2169 REAL ,DIMENSION(ims:ime,kms:kme,jms:jme,num_tracer),INTENT(INOUT) :: tracer_tend
2170 REAL ,DIMENSION(ims:ime,kms:kme,jms:jme,num_scalar),INTENT(INOUT) :: scalar
2171 REAL ,DIMENSION(ims:ime,kms:kme,jms:jme,num_scalar),INTENT(INOUT) :: scalar_tend
2173 INTEGER :: i,k,j, im
2174 INTEGER :: itf,ktf,jtf,itsu,jtsv
2176 !-----------------------------------------------------------------------
2178 !<DESCRIPTION>
2179 !
2180 ! calculate_phy_tend couples the physics tendencies to the column mass (mu),
2181 ! because prognostic equations are in flux form, but physics tendencies are
2182 ! computed for uncoupled variables.
2183 !
2184 !</DESCRIPTION>
2186 itf=MIN(ite,ide-1)
2187 jtf=MIN(jte,jde-1)
2188 ktf=MIN(kte,kde-1)
2189 itsu=MAX(its,ids+1)
2190 jtsv=MAX(jts,jds+1)
2192 ! radiation
2194 IF (config_flags%ra_lw_physics .gt. 0 .or. config_flags%ra_sw_physics .gt. 0) THEN
2196 DO J=jts,jtf
2197 DO K=kts,ktf
2198 DO I=its,itf
2199 RTHRATEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RTHRATEN(I,K,J)
2200 ENDDO
2201 ENDDO
2202 ENDDO
2204 ENDIF
2206 ! cumulus
2208 IF (config_flags%cu_physics .gt. 0) THEN
2210 DO J=jts,jtf
2211 DO I=its,itf
2212 DO K=kts,ktf
2213 RUCUTEN(I,K,J) =(c1(k)*mut(I,J)+c2(k))*RUCUTEN(I,K,J)
2214 RVCUTEN(I,K,J) =(c1(k)*mut(I,J)+c2(k))*RVCUTEN(I,K,J)
2215 RTHCUTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RTHCUTEN(I,K,J)
2216 RQVCUTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RQVCUTEN(I,K,J)
2217 ENDDO
2218 ENDDO
2219 ENDDO
2221 IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
2222 DO J=jts,jtf
2223 DO I=its,itf
2224 DO K=kts,ktf
2225 RQCCUTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RQCCUTEN(I,K,J)
2226 ENDDO
2227 ENDDO
2228 ENDDO
2229 ENDIF
2231 IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
2232 DO J=jts,jtf
2233 DO I=its,itf
2234 DO K=kts,ktf
2235 RQRCUTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RQRCUTEN(I,K,J)
2236 ENDDO
2237 ENDDO
2238 ENDDO
2239 ENDIF
2241 IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
2242 DO J=jts,jtf
2243 DO I=its,itf
2244 DO K=kts,ktf
2245 RQICUTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RQICUTEN(I,K,J)
2246 ENDDO
2247 ENDDO
2248 ENDDO
2249 ENDIF
2251 IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
2252 DO J=jts,jtf
2253 DO I=its,itf
2254 DO K=kts,ktf
2255 RQSCUTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RQSCUTEN(I,K,J)
2256 ENDDO
2257 ENDDO
2258 ENDDO
2259 ENDIF
2261 ENDIF
2263 ! shallow cumulus
2265 IF (config_flags%shcu_physics .gt. 0) THEN
2267 DO J=jts,jtf
2268 DO I=its,itf
2269 DO K=kts,ktf
2270 RUSHTEN(I,K,J) =(c1(k)*mut(I,J)+c2(k))*RUSHTEN(I,K,J)
2271 RVSHTEN(I,K,J) =(c1(k)*mut(I,J)+c2(k))*RVSHTEN(I,K,J)
2272 RTHSHTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RTHSHTEN(I,K,J)
2273 RQVSHTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RQVSHTEN(I,K,J)
2274 ENDDO
2275 ENDDO
2276 ENDDO
2278 IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
2279 DO J=jts,jtf
2280 DO I=its,itf
2281 DO K=kts,ktf
2282 RQCSHTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RQCSHTEN(I,K,J)
2283 ENDDO
2284 ENDDO
2285 ENDDO
2286 ENDIF
2288 IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
2289 DO J=jts,jtf
2290 DO I=its,itf
2291 DO K=kts,ktf
2292 RQRSHTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RQRSHTEN(I,K,J)
2293 ENDDO
2294 ENDDO
2295 ENDDO
2296 ENDIF
2298 IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
2299 DO J=jts,jtf
2300 DO I=its,itf
2301 DO K=kts,ktf
2302 RQISHTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RQISHTEN(I,K,J)
2303 ENDDO
2304 ENDDO
2305 ENDDO
2306 ENDIF
2308 IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
2309 DO J=jts,jtf
2310 DO I=its,itf
2311 DO K=kts,ktf
2312 RQSSHTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RQSSHTEN(I,K,J)
2313 ENDDO
2314 ENDDO
2315 ENDDO
2316 ENDIF
2318 IF(P_QG .ge. PARAM_FIRST_SCALAR)THEN
2319 DO J=jts,jtf
2320 DO I=its,itf
2321 DO K=kts,ktf
2322 RQGSHTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RQGSHTEN(I,K,J)
2323 ENDDO
2324 ENDDO
2325 ENDDO
2326 ENDIF
2328 ENDIF
2330 ! pbl
2332 IF (config_flags%bl_pbl_physics .gt. 0) THEN
2334 DO J=jts,jtf
2335 DO K=kts,ktf
2336 DO I=its,itf
2337 RUBLTEN(I,K,J) =(c1(k)*mut(I,J)+c2(k))*RUBLTEN(I,K,J)
2338 RVBLTEN(I,K,J) =(c1(k)*mut(I,J)+c2(k))*RVBLTEN(I,K,J)
2339 RTHBLTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RTHBLTEN(I,K,J)
2340 ENDDO
2341 ENDDO
2342 ENDDO
2344 IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
2345 DO J=jts,jtf
2346 DO K=kts,ktf
2347 DO I=its,itf
2348 RQVBLTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RQVBLTEN(I,K,J)
2349 ENDDO
2350 ENDDO
2351 ENDDO
2352 ENDIF
2354 IF (P_QC .ge. PARAM_FIRST_SCALAR) THEN
2355 DO J=jts,jtf
2356 DO K=kts,ktf
2357 DO I=its,itf
2358 RQCBLTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RQCBLTEN(I,K,J)
2359 ENDDO
2360 ENDDO
2361 ENDDO
2362 ENDIF
2364 IF (P_QI .ge. PARAM_FIRST_SCALAR) THEN
2365 DO J=jts,jtf
2366 DO K=kts,ktf
2367 DO I=its,itf
2368 RQIBLTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RQIBLTEN(I,K,J)
2369 ENDDO
2370 ENDDO
2371 ENDDO
2372 ENDIF
2374 ENDIF
2376 ! fdda
2377 ! note fdda u and v tendencies are staggered, also only interior points have muu/muv,
2378 ! so only couple those
2380 IF (config_flags%grid_fdda .gt. 0) THEN
2382 DO J=jts,jtf
2383 DO K=kts,ktf
2384 DO I=itsu,itf
2385 ! if( i == itf/2 .AND. j == jtf/2 .AND. k == ktf/2 ) &
2386 ! write(*,'(a,3i6,e15.5)') 'u_ten before=',i,k,j, RUNDGDTEN(i,k,j)
2387 RUNDGDTEN(I,K,J) =(c1(k)*muu(I,J)+c2(k))*RUNDGDTEN(I,K,J)
2388 ! if( i == itf/2 .AND. j == jtf/2 .AND. k==ktf/2 ) &
2389 ! write(*,'(a,2f15.5)') 'mu, muu=',(c1(k)*mut(I,J)+c2(k)), (c1(k)*muu(i,j)+c2(k))
2390 ! if( i == itf/2 .AND. j == jtf/2 .AND. k == ktf/2 ) &
2391 ! write(*,'(a,3i6,e15.5)') 'u_ten after=',i,k,j, RUNDGDTEN(i,k,j)
2392 ! if( RUNDGDTEN(i,k,j) > 30.0 ) write(*,*) 'IKJ=',i,k,j
2393 ENDDO
2394 ENDDO
2395 ENDDO
2396 ! write(*,'(a,e15.5)') 'u_ten MAXIMUM after=', maxval(RUNDGDTEN)
2397 DO J=jtsv,jtf
2398 DO K=kts,ktf
2399 DO I=its,itf
2400 RVNDGDTEN(I,K,J) =(c1(k)*muv(I,J)+c2(k))*RVNDGDTEN(I,K,J)
2401 ENDDO
2402 ENDDO
2403 ENDDO
2404 DO J=jts,jtf
2405 DO K=kts,ktf
2406 DO I=its,itf
2407 ! if( i == itf/2 .AND. j == jtf/2 .AND. k == ktf/2 ) &
2408 ! write(*,'(a,3i6,e15.5)') 'th before=',i,k,j, RTHNDGDTEN(I,K,J)
2409 RTHNDGDTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RTHNDGDTEN(I,K,J)
2410 ! RMUNDGDTEN(I,J) - no coupling
2411 ! if( i == itf/2 .AND. j == jtf/2 .AND. k == ktf/2 ) &
2412 ! write(*,'(a,3i6,e15.5)') 'th after=',i,k,j, RTHNDGDTEN(I,K,J)
2413 ENDDO
2414 ENDDO
2415 ENDDO
2417 IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
2418 DO J=jts,jtf
2419 DO K=kts,ktf
2420 DO I=its,itf
2421 RQVNDGDTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*RQVNDGDTEN(I,K,J)
2422 ENDDO
2423 ENDDO
2424 ENDDO
2425 ENDIF
2427 ENDIF
2429 ! 4d couple scalar tendencies that have only uncoupled physics tendencies at this point
2431 DO im = PARAM_FIRST_SCALAR,num_scalar
2432 DO J=jts,jtf
2433 DO K=kts,ktf
2434 DO I=its,itf
2435 scalar_tend(I,K,J,im)=(c1(k)*mut(I,J)+c2(k))*scalar_tend(I,K,J,im)
2436 ENDDO
2437 ENDDO
2438 ENDDO
2439 ENDDO
2441 DO im = PARAM_FIRST_SCALAR,num_tracer
2442 DO J=jts,jtf
2443 DO K=kts,ktf
2444 DO I=its,itf
2445 tracer_tend(I,K,J,im)=(c1(k)*mut(I,J)+c2(k))*tracer_tend(I,K,J,im)
2446 ENDDO
2447 ENDDO
2448 ENDDO
2449 ENDDO
2452 END SUBROUTINE calculate_phy_tend
2454 !-----------------------------------------------------------------------
2456 SUBROUTINE positive_definite_filter ( a, &
2457 ids,ide, jds,jde, kds,kde, &
2458 ims,ime, jms,jme, kms,kme, &
2459 its,ite, jts,jte, kts,kte )
2461 IMPLICIT NONE
2463 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
2464 ims,ime, jms,jme, kms,kme, &
2465 its,ite, jts,jte, kts,kte
2467 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: a
2469 INTEGER :: i,k,j
2471 !<DESCRIPTION>
2472 !
2473 ! debug and testing code for bounding a variable
2474 !
2475 !</DESCRIPTION>
2477 DO j=jts,min(jte,jde-1)
2478 DO k=kts,kte-1
2479 DO i=its,min(ite,ide-1)
2480 ! a(i,k,j) = max(a(i,k,j),0.)
2481 a(i,k,j) = min(1000.,max(a(i,k,j),0.))
2482 ENDDO
2483 ENDDO
2484 ENDDO
2486 END SUBROUTINE positive_definite_filter
2488 !-----------------------------------------------------------------------
2490 SUBROUTINE bound_tke ( tke, tke_upper_bound, &
2491 ids,ide, jds,jde, kds,kde, &
2492 ims,ime, jms,jme, kms,kme, &
2493 its,ite, jts,jte, kts,kte )
2495 IMPLICIT NONE
2497 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
2498 ims,ime, jms,jme, kms,kme, &
2499 its,ite, jts,jte, kts,kte
2501 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: tke
2502 REAL, INTENT( IN) :: tke_upper_bound
2504 INTEGER :: i,k,j
2506 !<DESCRIPTION>
2507 !
2508 ! bounds tke between zero and tke_upper_bound.
2509 !
2510 !</DESCRIPTION>
2512 DO j=jts,min(jte,jde-1)
2513 DO k=kts,kte-1
2514 DO i=its,min(ite,ide-1)
2515 tke(i,k,j) = min(tke_upper_bound,max(tke(i,k,j),0.))
2516 ENDDO
2517 ENDDO
2518 ENDDO
2520 END SUBROUTINE bound_tke
2522 !-----------------------------------------------------------------------
2524 SUBROUTINE bound_qna ( qna, &
2525 ids,ide, jds,jde, kds,kde, &
2526 ims,ime, jms,jme, kms,kme, &
2527 its,ite, jts,jte, kts,kte )
2529 IMPLICIT NONE
2531 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
2532 ims,ime, jms,jme, kms,kme, &
2533 its,ite, jts,jte, kts,kte
2535 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: qna
2537 INTEGER :: i,k,j
2539 !<DESCRIPTION>
2540 !
2541 ! bounds aerosol between zero and current value
2542 !
2543 !</DESCRIPTION>
2545 DO j=jts,min(jte,jde-1)
2546 DO k=kts,kte-1
2547 DO i=its,min(ite,ide-1)
2548 qna(i,k,j) = max(qna(i,k,j),0.)
2549 ENDDO
2550 ENDDO
2551 ENDDO
2553 END SUBROUTINE bound_qna
2555 !----------------------------------------------------------------------------------
2556 !cyl: Implement the forward Lagrangian trajectory calculation in WRF
2557 !Chiaying Lee RSMAS/UM
2558 !----------------------------------------------------------------------------------
2559 subroutine trajectory ( grid,config_flags, &
2560 dt,itimestep,ru_m, rv_m, ww_m, &
2561 mut,muu,muv,c1h,c2h,c1f,c2f, &
2562 rdx, rdy, rdn, rdnw,rdzw, &
2563 traj_i,traj_j,traj_k, &
2564 traj_long,traj_lat, &
2565 xlong,xlat, &
2566 msft,msfu,msfv, &
2567 ids, ide, jds, jde, kds, kde, &
2568 ims, ime, jms, jme, kms, kme, &
2569 its, ite, jts, jte, kts, kte )
2571 !---------------------------------------------------------------------------------------------------
2573 use module_date_time
2574 use module_utility
2575 use module_domain, only : domain_clock_get
2576 use module_trajectory, only : traject, traj_cnt
2578 implicit none
2579 !---------------------------------------------------------------------------------------------------
2580 ! Subroutine trajectory calculates forward Lagrangian trajectory
2581 ! of the selecting points. The trajectory of each point is subjected to u, v, and w.
2582 ! (Lee and Chen 2013, MWR).
2583 !
2584 ! The trajectories is initialized with given longitude (degree), latitude (degree), and height(eta level).
2585 !
2586 !--traj_i grid number in x direction (on mass grid), float
2587 !--traj_j grid number in y direction (on mass grid), float
2588 !--traj_k grid number in z direction (on mass grid), float
2589 !--traj_long longitude of trajectories, float
2590 !--traj_lat longitude of trajectories, float
2591 !
2592 !--------------------------------------------------------------------------------------------------
2593 TYPE(proj_info) :: proj
2594 TYPE(domain), INTENT(IN) :: grid
2595 TYPE (grid_config_rec_type) , INTENT(IN) :: config_flags
2596 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2597 ims, ime, jms, jme, kms, kme, &
2598 its, ite, jts, jte, kts, kte
2599 INTEGER , INTENT(IN ) :: itimestep
2600 REAL , INTENT(IN ) :: rdx, &
2601 rdy, &
2602 dt
2603 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
2604 rdnw
2606 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: c1h, c2h, c1f, c2f
2608 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
2609 INTENT(IN ) :: ru_m, &
2610 rv_m, &
2611 ww_m
2612 REAL , DIMENSION( ims:ime , jms:jme ) , &
2613 INTENT(IN ) ::xlong, xlat
2615 real, dimension(1:config_flags%num_traj), intent(inout) :: traj_i,traj_j,traj_k
2616 real, dimension(1:config_flags%num_traj), intent(inout) :: traj_long,traj_lat
2617 real, dimension(ims:ime,kms:kme,jms:jme),intent(in) :: rdzw
2618 real, dimension(ims:ime,kms:kme,jms:jme)::u,v,w
2619 real, dimension(ims:ime,jms:jme),intent(in)::msft,msfu,msfv
2620 real, dimension(ims:ime,jms:jme),intent(in)::muu,muv,mut
2621 integer :: j,k
2622 integer :: i_beg,i_end
2623 integer :: i_traj,j_traj,k_traj,tjk
2624 integer :: dm
2625 real :: traj_u,traj_v,traj_w
2626 real :: rdx_grid,rdy_grid,rdz_grid
2627 real :: deltx, delty, deltz,ax
2628 real :: const1
2629 real :: temp_i,temp_j
2630 integer :: i_u,j_v,k_w
2631 real :: u_a,u_b,u_c,u_d,v_a,v_b,v_c,v_d,w_a,w_b,w_c,w_d
2632 real :: d_a,d_b,d_c,d_d
2633 real :: u_temp_upper,u_temp_lower
2634 real :: v_temp_upper,v_temp_lower
2635 real :: w_temp_upper,w_temp_lower
2636 real :: eta_old, eta_new
2637 integer :: keta, keta_temp
2638 character(len=19) :: current_timestr, next_timestr, wrk_timestr
2639 character(len=256) :: dbg_mes
2640 logical :: has_proj_map
2641 ! varalbe for map projectory
2642 real:: known_lat, known_lon
2643 TYPE (grid_config_rec_type) :: config_flags_temp
2644 TYPE(WRFU_Time) :: current_time, next_time
2645 TYPE(WRFU_Time) :: start_time, stop_time
2647 config_flags_temp = config_flags
2648 call trajmapproj(grid, config_flags_temp,proj)
2651 ! convert ru_m, rv_m and ww_m in u,v,w
2652 const1=1.0/2.0/sqrt(2.0)
2654 dm = grid%id
2656 write(dbg_mes,'(''trajectory('',i2.2,''): Entering trajectory'')') dm
2657 call wrf_debug(200,trim(dbg_mes) )
2659 i_beg = max( 1,its-1 )
2660 i_end = min( ite+2,ide-1 )
2661 do j=max(1,jts-1),min(jte+2,jde-1)
2662 do k=kms,kme-1
2663 u(its:i_end,k,j)=ru_m(its:i_end,k,j)/(c1h(k)*MUU(its:i_end,j)+c2h(k))*msfu(its:i_end,j)
2664 v(its:i_end,k,j)=rv_m(its:i_end,k,j)/(c1h(k)*MUV(its:i_end,j)+c2h(k))*msfv(its:i_end,j)
2665 enddo
2666 do k=kms,kme
2667 w(its:i_end,k,j)=ww_m(its:i_end,k,j)/(c1f(k)*MUT(its:i_end,j)+c2f(k))*msft(its:i_end,j)
2668 enddo
2669 enddo
2670 has_proj_map = &
2671 ( proj%code .EQ. PROJ_LC ) .OR. &
2672 ( proj%code .EQ. PROJ_PS_WGS84 ) .OR. &
2673 ( proj%code .EQ. PROJ_ALBERS_NAD83 ) .OR. &
2674 ( proj%code .EQ. PROJ_MERC ) .OR. &
2675 ( proj%code .EQ. PROJ_LATLON ) .OR. &
2676 ( proj%code .EQ. PROJ_CYL ) .OR. &
2677 ( proj%code .EQ. PROJ_CASSINI ) .OR. &
2678 ( proj%code .EQ. PROJ_GAUSS ) .OR. &
2679 ( proj%code .EQ. PROJ_ROTLL )
2681 traj_loop: &
2682 do tjk = 1,traj_cnt(dm)
2683 !do tjk = 1,config_flags%num_traj
2684 traj_is_active: &
2685 if (traj_i(tjk) .ne. -9999.0) then
2686 eta_old = 0.0
2687 eta_new = 0.0
2688 keta = 0
2689 keta_temp = 0
2690 if( has_proj_map ) then
2691 call latlon_to_ij (proj, &
2692 traj_lat(tjk),traj_long(tjk),traj_i(tjk),traj_j(tjk))
2693 end if
2695 i_traj = floor(traj_i(tjk)) ! find the lower_left_bottom corner for
2696 !trajectory
2697 j_traj = floor(traj_j(tjk)) !
2698 k_traj = floor(traj_k(tjk)) !
2699 traj_in_domain: &
2700 if ((i_traj .ge. its .and. i_traj .le. ite .and. i_traj .lt. ide) .and. &
2701 (j_traj .ge. jts .and. j_traj .le. jte .and. j_traj .lt. jde) .and. &
2702 (k_traj .le. kte .and. k_traj .lt. kde)) then
2703 !-----------------------------------------------------------------------------
2704 ! is trajectory in time interval?
2705 !-----------------------------------------------------------------------------
2706 call domain_clock_get( grid, current_time=current_time, &
2707 current_timestr=current_timestr)
2708 call geth_newdate( next_timestr, current_timestr, int(grid%dt) )
2709 call wrf_atotime( next_timestr, next_time )
2710 wrk_timestr(1:19) = traject(tjk,dm)%start_time(1:19)
2711 ! write(*,*) ' '
2712 ! write(*,*) traject(tjk,dm)
2713 ! write(dbg_mes,'(''trajectory('',i2,2,''): tjk,start_time = '',i3,1x,a)') &
2714 ! dm,tjk,wrk_timestr
2715 ! call wrf_debug( 200,trim(dbg_mes) )
2716 call wrf_atotime( wrk_timestr(1:19), start_time )
2717 wrk_timestr(1:19) = traject(tjk,dm)%stop_time(1:19)
2718 call wrf_atotime( wrk_timestr(1:19), stop_time )
2719 is_in_time_interval: &
2720 if( next_time .ge. start_time .and. next_time .le. stop_time &
2721 .and. .not. traject(tjk,dm)%is_stationary ) then
2722 ! for u : check x stagger
2723 if (traj_i(tjk)-real(floor(traj_i(tjk))) .ge. 0.5 ) then
2724 i_u=floor(traj_i(tjk)) + 1
2725 else
2726 i_u=floor(traj_i(tjk))
2727 endif
2728 ! for layer k_traj
2729 if (k_traj .ge. 1 ) then
2730 u_a=u(i_u ,k_traj,j_traj+1)
2731 u_b=u(i_u ,k_traj,j_traj )
2732 u_c=u(i_u+1,k_traj,j_traj )
2733 u_d=u(i_u+1,k_traj,j_traj+1)
2734 else
2735 u_a=0.0
2736 u_b=0.0
2737 u_c=0.0
2738 u_d=0.0
2739 endif
2740 d_a=abs((real(i_u+1)-(traj_i(tjk)+0.5))*(real(j_traj )-traj_j(tjk)))
2741 d_b=abs((real(i_u+1)-(traj_i(tjk)+0.5))*(real(j_traj+1)-traj_j(tjk)))
2742 d_c=abs((real(i_u )-(traj_i(tjk)+0.5))*(real(j_traj+1)-traj_j(tjk)))
2743 d_d=abs((real(i_u )-(traj_i(tjk)+0.5))*(real(j_traj )-traj_j(tjk)))
2744 u_temp_lower=(u_a*d_a+u_b*d_b+u_c*d_c+u_d*d_d)/(d_a+d_b+d_c+d_d)
2746 !for layer k_traj+1
2747 u_a=u(i_u ,k_traj+1,j_traj+1)
2748 u_b=u(i_u ,k_traj+1,j_traj )
2749 u_c=u(i_u+1,k_traj+1,j_traj )
2750 u_d=u(i_u+1,k_traj+1,j_traj+1)
2751 d_a=abs((real(i_u+1)-(traj_i(tjk)+0.5))*(real(j_traj )-traj_j(tjk)))
2752 d_b=abs((real(i_u+1)-(traj_i(tjk)+0.5))*(real(j_traj+1)-traj_j(tjk)))
2753 d_c=abs((real(i_u )-(traj_i(tjk)+0.5))*(real(j_traj+1)-traj_j(tjk)))
2754 d_d=abs((real(i_u )-(traj_i(tjk)+0.5))*(real(j_traj )-traj_j(tjk)))
2755 u_temp_upper=(u_a*d_a+u_b*d_b+u_c*d_c+u_d*d_d)/(d_a+d_b+d_c+d_d)
2757 traj_u=u_temp_upper*abs(real(k_traj)-traj_k(tjk))+u_temp_lower*abs(real(k_traj+1)-traj_k(tjk))
2759 ! for v: check y-stagger
2760 if (traj_j(tjk)-real(floor(traj_j(tjk))) .ge. 0.5 ) then
2761 j_v=floor(traj_j(tjk)) + 1
2762 else
2763 j_v=floor(traj_j(tjk))
2764 endif
2765 ! for layer k_traj
2766 if (k_traj .ge. 1 ) then
2767 v_a=v(i_traj ,k_traj,j_v+1)
2768 v_b=v(i_traj ,k_traj,j_v )
2769 v_c=v(i_traj+1,k_traj,j_v )
2770 v_d=v(i_traj+1,k_traj,j_v+1)
2771 else
2772 v_a=0.0
2773 v_b=0.0
2774 v_c=0.0
2775 v_d=0.0
2776 endif
2777 d_a=abs((real(i_traj+1)-traj_i(tjk))*(real(j_v )-(traj_j(tjk)+0.5)))
2778 d_b=abs((real(i_traj+1)-traj_i(tjk))*(real(j_v+1)-(traj_j(tjk)+0.5)))
2779 d_c=abs((real(i_traj )-traj_i(tjk))*(real(j_v+1)-(traj_j(tjk)+0.5)))
2780 d_d=abs((real(i_traj )-traj_i(tjk))*(real(j_v )-(traj_j(tjk)+0.5)))
2781 v_temp_lower=(v_a*d_a+v_b*d_b+v_c*d_c+v_d*d_d)/(d_a+d_b+d_c+d_d)
2782 !for layer k_traj+1
2783 v_a=v(i_traj ,k_traj+1,j_v+1)
2784 v_b=v(i_traj ,k_traj+1,j_v )
2785 v_c=v(i_traj+1,k_traj+1,j_v )
2786 v_d=v(i_traj+1,k_traj+1,j_v+1)
2787 d_a=abs((real(i_traj+1)-traj_i(tjk))*(real(j_v )-(traj_j(tjk)+0.5)))
2788 d_b=abs((real(i_traj+1)-traj_i(tjk))*(real(j_v+1)-(traj_j(tjk)+0.5)))
2789 d_c=abs((real(i_traj )-traj_i(tjk))*(real(j_v+1)-(traj_j(tjk)+0.5)))
2790 d_d=abs((real(i_traj )-traj_i(tjk))*(real(j_v )-(traj_j(tjk)+0.5)))
2791 v_temp_upper=(v_a*d_a+v_b*d_b+v_c*d_c+v_d*d_d)/(d_a+d_b+d_c+d_d)
2793 traj_v=v_temp_upper*abs(real(k_traj)-traj_k(tjk))+v_temp_lower*abs(real(k_traj+1)-traj_k(tjk))
2796 !for w: check for z-stagger
2797 if (traj_k(tjk)-real(floor(traj_k(tjk))) .ge. 0.5 ) then
2798 k_w=floor(traj_k(tjk)) + 1
2799 else
2800 k_w=floor(traj_k(tjk))
2801 endif
2802 !for layer j_traj
2803 if (k_w .ge. 1) then
2804 w_b=w(i_traj ,k_w ,j_traj)
2805 w_c=w(i_traj+1,k_w ,j_traj)
2806 else
2807 w_b=0.0
2808 w_c=0.0
2809 endif
2810 w_a=w(i_traj ,k_w+1,j_traj)
2811 w_d=w(i_traj+1,k_w+1,j_traj)
2812 d_a=abs((real(i_traj+1)-traj_i(tjk))*(real(k_w )-(traj_k(tjk)+0.5)))
2813 d_b=abs((real(i_traj+1)-traj_i(tjk))*(real(k_w+1)-(traj_k(tjk)+0.5)))
2814 d_c=abs((real(i_traj )-traj_i(tjk))*(real(k_w+1)-(traj_k(tjk)+0.5)))
2815 d_d=abs((real(i_traj )-traj_i(tjk))*(real(k_w )-(traj_k(tjk)+0.5)))
2816 w_temp_lower=(w_a*d_a+w_b*d_b+w_c*d_c+w_d*d_d)/(d_a+d_b+d_c+d_d)
2817 !for layer j_traj+1
2818 if (k_w .ge. 1) then
2819 w_b=w(i_traj ,k_w ,j_traj+1)
2820 w_c=w(i_traj+1,k_w ,j_traj+1)
2821 else
2822 w_b=0.0
2823 w_c=0.0
2824 endif
2825 w_a=w(i_traj ,k_w+1,j_traj+1)
2826 w_d=w(i_traj+1,k_w+1,j_traj+1)
2827 d_a=abs((real(i_traj+1)-traj_i(tjk))*(real(k_w )-(traj_k(tjk)+0.5)))
2828 d_b=abs((real(i_traj+1)-traj_i(tjk))*(real(k_w+1)-(traj_k(tjk)+0.5)))
2829 d_c=abs((real(i_traj )-traj_i(tjk))*(real(k_w+1)-(traj_k(tjk)+0.5)))
2830 d_d=abs((real(i_traj )-traj_i(tjk))*(real(k_w )-(traj_k(tjk)+0.5)))
2831 w_temp_upper=(w_a*d_a+w_b*d_b+w_c*d_c+w_d*d_d)/(d_a+d_b+d_c+d_d)
2833 traj_w=w_temp_upper*abs(real(j_traj)-traj_j(tjk))+w_temp_lower*abs(real(j_traj+1)-traj_j(tjk))
2834 !get old eta
2835 eta_old=grid%znw(k_w+1)*abs(traj_k(tjk)+0.5-real(k_w))+grid%znw(k_w)*abs(traj_k(tjk)+0.5-real(k_w+1))
2837 rdx_grid=rdx*msft(i_traj,j_traj)! 1/(dx/msft)
2838 rdy_grid=rdy*msft(i_traj,j_traj)
2839 rdz_grid=rdnw(k_traj)
2840 deltx=traj_u*DT
2841 delty=traj_v*DT
2842 deltz=traj_w*DT
2843 traj_i(tjk)=traj_i(tjk)+deltx*rdx_grid
2844 traj_j(tjk)=traj_j(tjk)+delty*rdy_grid
2845 eta_new=eta_old+deltz
2846 ! get new traj_k(tjk)
2847 keta_temp = 0
2848 do keta=1, kme-1
2849 if (eta_new .le. grid%znw(keta) .and. eta_new .gt. grid%znw(keta+1)) then
2850 keta_temp=keta
2851 endif
2852 enddo
2853 if (keta_temp .eq. 0) then
2854 traj_k(tjk) = traj_k(tjk)
2855 else
2856 traj_k(tjk) = (real(keta_temp)*abs(eta_new-grid%znw(keta_temp+1))+ &
2857 real(keta_temp+1)*abs(eta_new-grid%znw(keta_temp))) &
2858 /(grid%znw(keta_temp)-grid%znw(keta_temp+1))
2859 traj_k(tjk) = traj_k(tjk)-0.5
2860 endif
2861 !! convert i,j,k into lon, lat
2862 if( has_proj_map ) then
2863 call ij_to_latlon (proj, traj_i(tjk), &
2864 traj_j(tjk),traj_lat(tjk),traj_long(tjk))
2865 endif
2866 endif is_in_time_interval
2867 else traj_in_domain
2868 traj_i(tjk) = -9999.0
2869 traj_j(tjk) = -9999.0
2870 traj_k(tjk) = -9999.0
2871 traj_long(tjk) = -9999.0
2872 traj_lat(tjk) = -9999.0
2873 endif traj_in_domain
2874 endif traj_is_active
2875 traj_i(tjk) = wrf_dm_max_real(traj_i(tjk))
2876 traj_j(tjk) = wrf_dm_max_real(traj_j(tjk))
2877 traj_k(tjk) = wrf_dm_max_real(traj_k(tjk))
2878 traj_long(tjk) = wrf_dm_max_real(traj_long(tjk))
2879 traj_lat(tjk) = wrf_dm_max_real(traj_lat(tjk))
2880 enddo traj_loop
2881 !end trajectory
2883 END SUBROUTINE trajectory
2884 !------------------------------------------------------------------------
2885 !cyl:Implement the map prpjection code for trajectory calculation
2886 ! Chiaying Lee RSMAS/UM
2887 !------------------------------------------------------------------------
2888 subroutine trajmapproj (grid,config_flags,ts_proj)
2889 !------------------------------------------------------------------------
2890 !!! This code calculates map-projection of trajectories. It is from share/wrf_timeseries.F
2891 !------------------------------------------------------------------------
2893 IMPLICIT NONE
2896 ! Arguments
2897 TYPE (domain), INTENT(IN) :: grid
2898 TYPE (grid_config_rec_type) , INTENT(IN) :: config_flags
2899 ! Externals
2900 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
2901 INTEGER, EXTERNAL :: get_unused_unit
2904 ! Local variables
2905 INTEGER :: ntsloc_temp
2906 INTEGER :: i, k, iunit
2907 REAL :: ts_rx, ts_ry, ts_xlat, ts_xlong, ts_hgt
2908 REAL :: known_lat, known_lon
2909 CHARACTER (LEN=132) :: message
2910 TYPE (PROJ_INFO), INTENT(out) :: ts_proj
2913 INTEGER :: ids, ide, jds, jde, kds, kde, &
2914 ims, ime, jms, jme, kms, kme, &
2915 ips, ipe, jps, jpe, kps, kpe, &
2916 imsx, imex, jmsx, jmex, kmsx, kmex, &
2917 ipsx, ipex, jpsx, jpex, kpsx, kpex, &
2918 imsy, imey, jmsy, jmey, kmsy, kmey, &
2919 ipsy, ipey, jpsy, jpey, kpsy, kpey
2920 TYPE (grid_config_rec_type) :: config_flags_temp
2923 config_flags_temp = config_flags
2925 CALL get_ijk_from_grid ( grid , &
2926 ids, ide, jds, jde, kds, kde, &
2927 ims, ime, jms, jme, kms, kme, &
2928 ips, ipe, jps, jpe, kps, kpe, &
2929 imsx, imex, jmsx, jmex, kmsx, kmex, &
2930 ipsx, ipex, jpsx, jpex, kpsx, kpex, &
2931 imsy, imey, jmsy, jmey, kmsy, kmey, &
2932 ipsy, ipey, jpsy, jpey, kpsy, kpey )
2934 CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags_temp )
2937 ! Set up map transformation structure
2938 CALL map_init(ts_proj)
2941 IF (ips <= 1 .AND. 1 <= ipe .AND. &
2942 jps <= 1 .AND. 1 <= jpe) THEN
2943 known_lat = grid%xlat(1,1)
2944 known_lon = grid%xlong(1,1)
2945 ELSE
2946 known_lat = 9999.
2947 known_lon = 9999.
2948 END IF
2949 known_lat = wrf_dm_min_real(known_lat)
2950 known_lon = wrf_dm_min_real(known_lon)
2953 ! Mercator
2954 IF (config_flags%map_proj == PROJ_MERC) THEN
2955 CALL map_set(PROJ_MERC, ts_proj, &
2956 truelat1 = config_flags%truelat1, &
2957 lat1 = known_lat, &
2958 lon1 = known_lon, &
2959 knowni = 1., &
2960 knownj = 1., &
2961 dx = config_flags%dx)
2962 ! Lambert conformal
2963 ELSE IF (config_flags%map_proj == PROJ_LC) THEN
2964 CALL map_set(PROJ_LC, ts_proj, &
2965 truelat1 = config_flags%truelat1, &
2966 truelat2 = config_flags%truelat2, &
2967 stdlon = config_flags%stand_lon, &
2968 lat1 = known_lat, &
2969 lon1 = known_lon, &
2970 knowni = 1., &
2971 knownj = 1., &
2972 dx = config_flags%dx)
2973 ! Polar stereographic
2974 ELSE IF (config_flags%map_proj == PROJ_PS) THEN
2975 CALL map_set(PROJ_PS, ts_proj, &
2976 truelat1 = config_flags%truelat1, &
2977 stdlon = config_flags%stand_lon, &
2978 lat1 = known_lat, &
2979 lon1 = known_lon, &
2980 knowni = 1., &
2981 knownj = 1., &
2982 dx = config_flags%dx)
2985 ! Cassini (global ARW)
2986 ELSE IF (config_flags%map_proj == PROJ_CASSINI) THEN
2987 CALL map_set(PROJ_CASSINI, ts_proj, &
2988 latinc = grid%dy*360.0/(2.0*EARTH_RADIUS_M*PI), &
2989 loninc = grid%dx*360.0/(2.0*EARTH_RADIUS_M*PI), &
2990 lat1 = known_lat, &
2991 lon1 = known_lon, &
2992 ! We still need to get POLE_LAT and POLE_LON metadata variables before
2993 ! this will work for rotated poles.
2994 lat0 = 90.0, &
2995 lon0 = 0.0, &
2996 knowni = 1., &
2997 knownj = 1., &
2998 stdlon = config_flags%stand_lon)
3001 ! Rotated latitude-longitude
3002 ELSE IF (config_flags%map_proj == PROJ_ROTLL) THEN
3003 CALL map_set(PROJ_ROTLL, ts_proj, &
3004 ! I have no idea how this should work for NMM nested domains
3005 ixdim = grid%e_we-1, &
3006 jydim = grid%e_sn-1, &
3007 phi = real(grid%e_sn-2)*grid%dy/2.0, &
3008 lambda = real(grid%e_we-2)*grid%dx, &
3009 lat1 = config_flags%cen_lat, &
3010 lon1 = config_flags%cen_lon, &
3011 latinc = grid%dy, &
3012 loninc = grid%dx, &
3013 stagger = HH)
3016 END IF
3019 end subroutine trajmapproj
3021 END MODULE module_em