Context Navigation

module_em.F @ 3553

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

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

Download in other formats:

Original Format