Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_em.F @ 1393

Last change on this file since 1393 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 77.8 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format