Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_em.F @ 3567

Last change on this file since 3567 was 2761, checked in by aslmd, 2 years ago
Applied planetary adaptation changes to WRFV3. job done previously by LMD_LES_MARS_install. Moved Registry.EM.newphys to put it simply in Registry.EM
File size: 71.4 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	xlat, f, e, sina, cosa, &
182	fnm, fnp, rdn, rdnw, &
183	dt, rdx, rdy, khdif, kvdif, xkmhd, xkhh, &
184	diff_6th_opt, diff_6th_factor, &
185	dampcoef,zdamp,damp_opt, &
186	cf1, cf2, cf3, cfn, cfn1, n_moist, &
187	non_hydrostatic, 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	xlat, &
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
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) &
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
474	CALL w_damp ( rw_tend, max_vert_cfl, &
475	max_horiz_cfl, &
476	u, v, ww, w, mut, rdnw, &
477	rdx, rdy, msfux, msfuy, msfvx, &
478	msfvy, dt, config_flags, &
479	ids, ide, jds, jde, kds, kde, &
480	ims, ime, jms, jme, kms, kme, &
481	its, ite, jts, jte, kts, kte )
482
483	IF(config_flags%pert_coriolis) THEN
484
485	CALL perturbation_coriolis ( ru, rv, rw, &
486	ru_tend, rv_tend, rw_tend, &
487	config_flags, &
488	u_base, v_base, z_base, &
489	muu, muv, phb, ph, &
490	msftx, msfty, msfux, msfuy, &
491	msfvx, msfvy, &
492	f, e, sina, cosa, fnm, fnp, &
493	ids, ide, jds, jde, kds, kde, &
494	ims, ime, jms, jme, kms, kme, &
495	its, ite, jts, jte, kts, kte )
496	ELSE
497	CALL coriolis ( ru, rv, rw, &
498	ru_tend, rv_tend, rw_tend, &
499	config_flags, &
500	msftx, msfty, msfux, msfuy, &
501	msfvx, msfvy, &
502	f, e, sina, cosa, fnm, fnp, &
503	ids, ide, jds, jde, kds, kde, &
504	ims, ime, jms, jme, kms, kme, &
505	its, ite, jts, jte, kts, kte )
506
507	END IF
508
509	CALL curvature ( ru, rv, rw, u, v, w, &
510	ru_tend, rv_tend, rw_tend, &
511	config_flags, &
512	msfux, msfuy, msfvx, msfvy, &
513	msftx, msfty, &
514	xlat, fnm, fnp, rdx, rdy, &
515	ids, ide, jds, jde, kds, kde, &
516	ims, ime, jms, jme, kms, kme, &
517	its, ite, jts, jte, kts, kte )
518
519	! Damping option added for Held-Suarez test (also uses lw option HELDSUAREZ)
520	IF (config_flags%ra_lw_physics == HELDSUAREZ) THEN
521	CALL held_suarez_damp ( ru_tend, rv_tend, &
522	ru,rv,p,pb, &
523	ids, ide, jds, jde, kds, kde, &
524	ims, ime, jms, jme, kms, kme, &
525	its, ite, jts, jte, kts, kte )
526	END IF
527
528	!**************************************************************
529	!
530	! Next, the terms that we integrate only with forward-in-time
531	! (evaluate with time t variables).
532	!
533	!**************************************************************
534
535	forward_step: IF( rk_step == 1 ) THEN
536
537	diff_opt1 : IF (config_flags%diff_opt .eq. 1) THEN
538
539	CALL horizontal_diffusion ('u', u, ru_tendf, mut, config_flags, &
540	msfux, msfuy, msfvx, msfvx_inv, &
541	msfvy,msftx, msfty, &
542	khdif, xkmhd, rdx, rdy, &
543	ids, ide, jds, jde, kds, kde, &
544	ims, ime, jms, jme, kms, kme, &
545	its, ite, jts, jte, kts, kte )
546
547	CALL horizontal_diffusion ('v', v, rv_tendf, mut, config_flags, &
548	msfux, msfuy, msfvx, msfvx_inv, &
549	msfvy,msftx, msfty, &
550	khdif, xkmhd, rdx, rdy, &
551	ids, ide, jds, jde, kds, kde, &
552	ims, ime, jms, jme, kms, kme, &
553	its, ite, jts, jte, kts, kte )
554
555	CALL horizontal_diffusion ('w', w, rw_tendf, mut, config_flags, &
556	msfux, msfuy, msfvx, msfvx_inv, &
557	msfvy,msftx, msfty, &
558	khdif, xkmhd, rdx, rdy, &
559	ids, ide, jds, jde, kds, kde, &
560	ims, ime, jms, jme, kms, kme, &
561	its, ite, jts, jte, kts, kte )
562
563	khdq = 3.*khdif
564	CALL horizontal_diffusion_3dmp ( 'm', t, t_tendf, mut, &
565	config_flags, t_init, &
566	msfux, msfuy, msfvx, msfvx_inv, &
567	msfvy, msftx, msfty, &
568	khdq , xkhh, rdx, rdy, &
569	ids, ide, jds, jde, kds, kde, &
570	ims, ime, jms, jme, kms, kme, &
571	its, ite, jts, jte, kts, kte )
572
573	! pbl_test : IF (config_flags%bl_pbl_physics .eq. 0) THEN
574	!!!****MARS: vertical diffusion is done in the physics (TODO: consider the nonhydrostatic case ?)
575	pbl_test : IF ( (config_flags%bl_pbl_physics .eq. 0) &
576	.AND. (.not. config_flags%modif_wrf) ) THEN
577
578	CALL vertical_diffusion_u ( u, ru_tendf, config_flags, &
579	u_base, &
580	alt, muu, rdn, rdnw, kvdif, &
581	ids, ide, jds, jde, kds, kde, &
582	ims, ime, jms, jme, kms, kme, &
583	its, ite, jts, jte, kts, kte )
584
585	CALL vertical_diffusion_v ( v, rv_tendf, config_flags, &
586	v_base, &
587	alt, muv, rdn, rdnw, kvdif, &
588	ids, ide, jds, jde, kds, kde, &
589	ims, ime, jms, jme, kms, kme, &
590	its, ite, jts, jte, kts, kte )
591
592	IF (non_hydrostatic) &
593	CALL vertical_diffusion ( 'w', w, rw_tendf, config_flags, &
594	alt, mut, rdn, rdnw, kvdif, &
595	ids, ide, jds, jde, kds, kde, &
596	ims, ime, jms, jme, kms, kme, &
597	its, ite, jts, jte, kts, kte )
598
599	kvdq = 3.*kvdif
600	CALL vertical_diffusion_3dmp ( t, t_tendf, config_flags, t_init, &
601	alt, mut, rdn, rdnw, kvdq , &
602	ids, ide, jds, jde, kds, kde, &
603	ims, ime, jms, jme, kms, kme, &
604	its, ite, jts, jte, kts, kte )
605
606	ENDIF pbl_test
607
608	! Theta tendency computations.
609
610	END IF diff_opt1
611
612	IF ( diff_6th_opt .NE. 0 ) THEN
613
614	CALL sixth_order_diffusion( 'u', u, ru_tendf, mut, dt, &
615	config_flags, &
616	diff_6th_opt, diff_6th_factor, &
617	ids, ide, jds, jde, kds, kde, &
618	ims, ime, jms, jme, kms, kme, &
619	its, ite, jts, jte, kts, kte )
620
621	CALL sixth_order_diffusion( 'v', v, rv_tendf, mut, dt, &
622	config_flags, &
623	diff_6th_opt, diff_6th_factor, &
624	ids, ide, jds, jde, kds, kde, &
625	ims, ime, jms, jme, kms, kme, &
626	its, ite, jts, jte, kts, kte )
627
628	IF (non_hydrostatic) &
629	CALL sixth_order_diffusion( 'w', w, rw_tendf, mut, dt, &
630	config_flags, &
631	diff_6th_opt, diff_6th_factor, &
632	ids, ide, jds, jde, kds, kde, &
633	ims, ime, jms, jme, kms, kme, &
634	its, ite, jts, jte, kts, kte )
635
636	CALL sixth_order_diffusion( 'm', t, t_tendf, mut, dt, &
637	config_flags, &
638	diff_6th_opt, diff_6th_factor, &
639	ids, ide, jds, jde, kds, kde, &
640	ims, ime, jms, jme, kms, kme, &
641	its, ite, jts, jte, kts, kte )
642
643	ENDIF
644
645	IF( damp_opt .eq. 2 ) &
646	CALL rk_rayleigh_damp( ru_tendf, rv_tendf, &
647	rw_tendf, t_tendf, &
648	u, v, w, t, t_init, &
649	mut, muu, muv, ph, phb, &
650	u_base, v_base, t_base, z_base, &
651	dampcoef, zdamp, &
652	ids, ide, jds, jde, kds, kde, &
653	ims, ime, jms, jme, kms, kme, &
654	its, ite, jts, jte, kts, kte )
655
656	END IF forward_step
657
658	END SUBROUTINE rk_tendency
659
660	!-------------------------------------------------------------------------------
661
662	SUBROUTINE rk_addtend_dry ( ru_tend, rv_tend, rw_tend, ph_tend, t_tend, &
663	ru_tendf, rv_tendf, rw_tendf, ph_tendf, t_tendf, &
664	u_save, v_save, w_save, ph_save, t_save, &
665	mu_tend, mu_tendf, rk_step, &
666	h_diabatic, mut, msftx, msfty, msfux, msfuy, &
667	msfvx, msfvx_inv, msfvy, &
668	ids,ide, jds,jde, kds,kde, &
669	ims,ime, jms,jme, kms,kme, &
670	ips,ipe, jps,jpe, kps,kpe, &
671	its,ite, jts,jte, kts,kte )
672
673	IMPLICIT NONE
674
675	! Input data.
676
677	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
678	ims, ime, jms, jme, kms, kme, &
679	ips, ipe, jps, jpe, kps, kpe, &
680	its, ite, jts, jte, kts, kte
681	INTEGER , INTENT(IN ) :: rk_step
682
683	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(INOUT) :: ru_tend, &
684	rv_tend, &
685	rw_tend, &
686	ph_tend, &
687	t_tend, &
688	ru_tendf, &
689	rv_tendf, &
690	rw_tendf, &
691	ph_tendf, &
692	t_tendf
693
694	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(INOUT) :: mu_tend, &
695	mu_tendf
696
697	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(IN ) :: u_save, &
698	v_save, &
699	w_save, &
700	ph_save, &
701	t_save, &
702	h_diabatic
703
704	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut, &
705	msftx, &
706	msfty, &
707	msfux, &
708	msfuy, &
709	msfvx, &
710	msfvx_inv, &
711	msfvy
712
713
714	! Local
715	INTEGER :: i, j, k
716
717	!<DESCRIPTION>
718	!
719	! rk_addtend_dry constructs the full large-timestep tendency terms for
720	! momentum (u,v,w), theta and geopotential equations. This is accomplished
721	! by combining the physics tendencies (in *tendf; these are computed
722	! the first RK substep, held fixed thereafter) with the RK tendencies
723	! (in *tend, these include advection, pressure gradient, etc;
724	! these change each rk substep). Output is in *tend.
725	!
726	!</DESCRIPTION>
727
728	! Finally, add the forward-step tendency to the rk_tendency
729
730	! u/v/w/save contain bc tendency that needs to be multiplied by msf
731	! (u by msfuy, v by msfvx)
732	! before adding it to physics tendency (*tendf)
733	! For momentum we need the final tendency to include an inverse msf
734	! physics/bc tendency needs to be divided, advection tendency already has it
735
736	! For scalars we need the final tendency to include an inverse msf (msfty)
737	! advection tendency is OK, physics/bc tendency needs to be divided by msf
738
739	DO j = jts,MIN(jte,jde-1)
740	DO k = kts,kte-1
741	DO i = its,ite
742	! multiply by my to uncouple u
743	IF(rk_step == 1)ru_tendf(i,k,j) = ru_tendf(i,k,j) + u_save(i,k,j)*msfuy(i,j)
744	! divide by my to couple u
745	ru_tend(i,k,j) = ru_tend(i,k,j) + ru_tendf(i,k,j)/msfuy(i,j)
746	ENDDO
747	ENDDO
748	ENDDO
749
750	DO j = jts,jte
751	DO k = kts,kte-1
752	DO i = its,MIN(ite,ide-1)
753	! multiply by mx to uncouple v
754	IF(rk_step == 1)rv_tendf(i,k,j) = rv_tendf(i,k,j) + v_save(i,k,j)*msfvx(i,j)
755	! divide by mx to couple v
756	rv_tend(i,k,j) = rv_tend(i,k,j) + rv_tendf(i,k,j)*msfvx_inv(i,j)
757	ENDDO
758	ENDDO
759	ENDDO
760
761	DO j = jts,MIN(jte,jde-1)
762	DO k = kts,kte
763	DO i = its,MIN(ite,ide-1)
764	! multiply by my to uncouple w
765	IF(rk_step == 1)rw_tendf(i,k,j) = rw_tendf(i,k,j) + w_save(i,k,j)*msfty(i,j)
766	! divide by my to couple w
767	rw_tend(i,k,j) = rw_tend(i,k,j) + rw_tendf(i,k,j)/msfty(i,j)
768	IF(rk_step == 1)ph_tendf(i,k,j) = ph_tendf(i,k,j) + ph_save(i,k,j)
769	! divide by my to couple scalar
770	ph_tend(i,k,j) = ph_tend(i,k,j) + ph_tendf(i,k,j)/msfty(i,j)
771	ENDDO
772	ENDDO
773	ENDDO
774
775	DO j = jts,MIN(jte,jde-1)
776	DO k = kts,kte-1
777	DO i = its,MIN(ite,ide-1)
778	IF(rk_step == 1)t_tendf(i,k,j) = t_tendf(i,k,j) + t_save(i,k,j)
779	! divide by my to couple theta
780	t_tend(i,k,j) = t_tend(i,k,j) + t_tendf(i,k,j)/msfty(i,j) &
781	+ mut(i,j)*h_diabatic(i,k,j)/msfty(i,j)
782	! divide by my to couple heating
783	ENDDO
784	ENDDO
785	ENDDO
786
787	DO j = jts,MIN(jte,jde-1)
788	DO i = its,MIN(ite,ide-1)
789	! mu tendencies not coupled with 1/msf
790	mu_tend(i,j) = mu_tend(i,j) + mu_tendf(i,j)
791	ENDDO
792	ENDDO
793
794	END SUBROUTINE rk_addtend_dry
795
796	!-------------------------------------------------------------------------------
797
798	SUBROUTINE rk_scalar_tend ( scs, sce, config_flags, &
799	rk_step, dt, &
800	ru, rv, ww, mut, mub, mu_old, &
801	alt, &
802	scalar_old, scalar, &
803	scalar_tends, advect_tend, &
804	RQVFTEN, &
805	base, moist_step, fnm, fnp, &
806	msfux, msfuy, msfvx, msfvx_inv, &
807	msfvy, msftx, msfty, &
808	rdx, rdy, rdn, rdnw, &
809	khdif, kvdif, xkmhd, &
810	diff_6th_opt, diff_6th_factor, &
811	pd_advection, &
812	ids, ide, jds, jde, kds, kde, &
813	ims, ime, jms, jme, kms, kme, &
814	its, ite, jts, jte, kts, kte )
815
816	IMPLICIT NONE
817
818	! Input data.
819
820	TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
821
822	INTEGER , INTENT(IN ) :: rk_step, scs, sce
823	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
824	ims, ime, jms, jme, kms, kme, &
825	its, ite, jts, jte, kts, kte
826
827	LOGICAL , INTENT(IN ) :: moist_step
828
829	REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce ), &
830	INTENT(IN ) :: scalar, scalar_old
831
832	REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce ), &
833	INTENT(INOUT) :: scalar_tends
834
835	REAL, DIMENSION(ims:ime, kms:kme, jms:jme ), INTENT(INOUT) :: advect_tend
836
837	REAL, DIMENSION(ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RQVFTEN
838
839	REAL, DIMENSION(ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: ru, &
840	rv, &
841	ww, &
842	xkmhd, &
843	alt
844
845
846	REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm, &
847	fnp, &
848	rdn, &
849	rdnw, &
850	base
851
852	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
853	msfuy, &
854	msfvx, &
855	msfvx_inv, &
856	msfvy, &
857	msftx, &
858	msfty, &
859	mub, &
860	mut, &
861	mu_old
862
863	REAL , INTENT(IN ) :: rdx, &
864	rdy, &
865	khdif, &
866	kvdif
867
868	INTEGER, INTENT( IN ) :: diff_6th_opt
869	REAL, INTENT( IN ) :: diff_6th_factor
870
871	REAL , INTENT(IN ) :: dt
872
873	LOGICAL, INTENT(IN ) :: pd_advection
874
875	! Local data
876
877	INTEGER :: im, i,j,k
878	INTEGER :: time_step
879
880	REAL :: khdq, kvdq, tendency
881
882	!<DESCRIPTION>
883	!
884	! rk_scalar_tend calls routines that computes scalar tendency from advection
885	! and 3D mixing (TKE or fixed eddy viscosities).
886	!
887	!</DESCRIPTION>
888
889
890	khdq = khdif/prandtl
891	kvdq = kvdif/prandtl
892
893	scalar_loop : DO im = scs, sce
894
895	CALL zero_tend ( advect_tend(ims,kms,jms), &
896	ids, ide, jds, jde, kds, kde, &
897	ims, ime, jms, jme, kms, kme, &
898	its, ite, jts, jte, kts, kte )
899
900	CALL nl_get_time_step ( 1, time_step )
901
902	IF( (rk_step == 3) .and. pd_advection ) THEN
903
904	CALL advect_scalar_pd ( scalar(ims,kms,jms,im), &
905	scalar_old(ims,kms,jms,im), &
906	advect_tend(ims,kms,jms), &
907	ru, rv, ww, mut, mub, mu_old, &
908	config_flags, &
909	msfux, msfuy, msfvx, msfvy, &
910	msftx, msfty, fnm, fnp, &
911	rdx, rdy, rdnw,dt, &
912	ids, ide, jds, jde, kds, kde, &
913	ims, ime, jms, jme, kms, kme, &
914	its, ite, jts, jte, kts, kte )
915
916	ELSE
917
918	CALL advect_scalar ( scalar(ims,kms,jms,im), &
919	scalar(ims,kms,jms,im), &
920	advect_tend(ims,kms,jms), &
921	ru, rv, ww, mut, time_step, &
922	config_flags, &
923	msfux, msfuy, msfvx, msfvy, &
924	msftx, msfty, fnm, fnp, &
925	rdx, rdy, rdnw, &
926	ids, ide, jds, jde, kds, kde, &
927	ims, ime, jms, jme, kms, kme, &
928	its, ite, jts, jte, kts, kte )
929	END IF
930
931	IF((config_flags%cu_physics == GDSCHEME .OR. config_flags%cu_physics == G3SCHEME) &
932	.and. moist_step .and. ( im == P_QV) ) THEN
933
934	CALL set_tend( RQVFTEN, advect_tend, msfty, &
935	ids, ide, jds, jde, kds, kde, &
936	ims, ime, jms, jme, kms, kme, &
937	its, ite, jts, jte, kts, kte )
938	ENDIF
939
940	rk_step_1: IF( rk_step == 1 ) THEN
941
942	diff_opt1 : IF (config_flags%diff_opt .eq. 1) THEN
943
944	CALL horizontal_diffusion ( 'm', scalar(ims,kms,jms,im), &
945	scalar_tends(ims,kms,jms,im), mut, &
946	config_flags, &
947	msfux, msfuy, msfvx, msfvx_inv, &
948	msfvy, msftx, msfty, &
949	khdq , xkmhd, rdx, rdy, &
950	ids, ide, jds, jde, kds, kde, &
951	ims, ime, jms, jme, kms, kme, &
952	its, ite, jts, jte, kts, kte )
953
954	!!!****MARS: done in the physics
955	pbl_test : IF ( (config_flags%bl_pbl_physics .eq. 0) &
956	.AND. (.not. config_flags%modif_wrf) ) THEN
957	! pbl_test : IF (config_flags%bl_pbl_physics .eq. 0) THEN
958
959	IF( (moist_step) .and. ( im == P_QV)) THEN
960
961	CALL vertical_diffusion_mp ( scalar(ims,kms,jms,im), &
962	scalar_tends(ims,kms,jms,im), &
963	config_flags, base, &
964	alt, mut, rdn, rdnw, kvdq , &
965	ids, ide, jds, jde, kds, kde, &
966	ims, ime, jms, jme, kms, kme, &
967	its, ite, jts, jte, kts, kte )
968
969	ELSE
970
971	CALL vertical_diffusion ( 'm', scalar(ims,kms,jms,im), &
972	scalar_tends(ims,kms,jms,im), &
973	config_flags, &
974	alt, mut, rdn, rdnw, kvdq, &
975	ids, ide, jds, jde, kds, kde, &
976	ims, ime, jms, jme, kms, kme, &
977	its, ite, jts, jte, kts, kte )
978
979	END IF
980
981	ENDIF pbl_test
982
983	ENDIF diff_opt1
984
985	IF ( diff_6th_opt .NE. 0 ) &
986	CALL sixth_order_diffusion( 'm', scalar(ims,kms,jms,im), &
987	scalar_tends(ims,kms,jms,im), &
988	mut, dt, config_flags, &
989	diff_6th_opt, diff_6th_factor, &
990	ids, ide, jds, jde, kds, kde, &
991	ims, ime, jms, jme, kms, kme, &
992	its, ite, jts, jte, kts, kte )
993
994	ENDIF rk_step_1
995
996	END DO scalar_loop
997
998	END SUBROUTINE rk_scalar_tend
999
1000	!-------------------------------------------------------------------------------
1001
1002	SUBROUTINE rk_update_scalar( scs, sce, &
1003	scalar_1, scalar_2, sc_tend, &
1004	advect_tend, msftx, msfty, &
1005	mu_old, mu_new, mu_base, &
1006	rk_step, dt, spec_zone, &
1007	config_flags, &
1008	ids, ide, jds, jde, kds, kde, &
1009	ims, ime, jms, jme, kms, kme, &
1010	its, ite, jts, jte, kts, kte )
1011
1012	IMPLICIT NONE
1013
1014	! Input data.
1015
1016	TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
1017
1018	INTEGER , INTENT(IN ) :: scs, sce, rk_step, spec_zone
1019	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1020	ims, ime, jms, jme, kms, kme, &
1021	its, ite, jts, jte, kts, kte
1022
1023	REAL, INTENT(IN ) :: dt
1024
1025	REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce), &
1026	INTENT(INOUT) :: scalar_1, &
1027	scalar_2, &
1028	sc_tend
1029
1030	REAL, DIMENSION(ims:ime, kms:kme, jms:jme ), &
1031	INTENT(IN) :: advect_tend
1032
1033	REAL, DIMENSION(ims:ime, jms:jme ), INTENT(IN ) :: mu_old, &
1034	mu_new, &
1035	mu_base, &
1036	msftx, &
1037	msfty
1038
1039	INTEGER :: i,j,k,im
1040	REAL :: sc_middle, msfsq
1041	REAL, DIMENSION(its:ite) :: muold, r_munew
1042
1043	REAL, DIMENSION(its:ite, kts:kte, jts:jte ) :: tendency
1044
1045	INTEGER :: i_start,i_end,j_start,j_end,k_start,k_end
1046	INTEGER :: i_start_spc,i_end_spc,j_start_spc,j_end_spc,k_start_spc,k_end_spc
1047
1048	!<DESCRIPTION>
1049	!
1050	! rk_scalar_update advances the scalar equation given the time t value
1051	! of the scalar and the scalar tendency.
1052	!
1053	!</DESCRIPTION>
1054
1055
1056	!
1057	! set loop limits.
1058
1059	i_start = its
1060	i_end = ite
1061	j_start = jts
1062	j_end = jte
1063	k_start = kts
1064	k_end = kte-1
1065	IF(j_end == jde) j_end = j_end - 1
1066	IF(i_end == ide) i_end = i_end - 1
1067
1068	i_start_spc = i_start
1069	i_end_spc = i_end
1070	j_start_spc = j_start
1071	j_end_spc = j_end
1072	k_start_spc = k_start
1073	k_end_spc = k_end
1074
1075	IF( config_flags%nested .or. config_flags%specified ) THEN
1076	IF( .NOT. config_flags%periodic_x)i_start = max( its,ids+spec_zone )
1077	IF( .NOT. config_flags%periodic_x)i_end = min( ite,ide-spec_zone-1 )
1078	j_start = max( jts,jds+spec_zone )
1079	j_end = min( jte,jde-spec_zone-1 )
1080	k_start = kts
1081	k_end = min( kte, kde-1 )
1082	ENDIF
1083
1084	IF ( rk_step == 1 ) THEN
1085
1086	! replace t-dt values (in scalar_1) with t values scalar_2,
1087	! then compute new values by adding tendency to values at t
1088
1089	DO im = scs,sce
1090
1091	DO j = jts, min(jte,jde-1)
1092	DO k = kts, min(kte,kde-1)
1093	DO i = its, min(ite,ide-1)
1094	tendency(i,k,j) = 0.
1095	ENDDO
1096	ENDDO
1097	ENDDO
1098
1099	DO j = j_start,j_end
1100	DO k = k_start,k_end
1101	DO i = i_start,i_end
1102	! scalar was coupled with my
1103	tendency(i,k,j) = advect_tend(i,k,j) * msfty(i,j)
1104	ENDDO
1105	ENDDO
1106	ENDDO
1107
1108	DO j = j_start_spc,j_end_spc
1109	DO k = k_start_spc,k_end_spc
1110	DO i = i_start_spc,i_end_spc
1111	tendency(i,k,j) = tendency(i,k,j) + sc_tend(i,k,j,im)
1112	ENDDO
1113	ENDDO
1114	ENDDO
1115
1116	DO j = jts, min(jte,jde-1)
1117
1118	DO i = its, min(ite,ide-1)
1119	muold(i) = mu_old(i,j) + mu_base(i,j)
1120	r_munew(i) = 1./(mu_new(i,j) + mu_base(i,j))
1121	ENDDO
1122
1123	DO k = kts, min(kte,kde-1)
1124	DO i = its, min(ite,ide-1)
1125
1126	scalar_1(i,k,j,im) = scalar_2(i,k,j,im)
1127	scalar_2(i,k,j,im) = (muold(i)*scalar_1(i,k,j,im) &
1128	+ dttendency(i,k,j))r_munew(i)
1129
1130	ENDDO
1131	ENDDO
1132	ENDDO
1133
1134	ENDDO
1135
1136	ELSE
1137
1138	! just compute new values, scalar_1 already at time t.
1139
1140	DO im = scs, sce
1141
1142	DO j = jts, min(jte,jde-1)
1143	DO k = kts, min(kte,kde-1)
1144	DO i = its, min(ite,ide-1)
1145	tendency(i,k,j) = 0.
1146	ENDDO
1147	ENDDO
1148	ENDDO
1149
1150	DO j = j_start,j_end
1151	DO k = k_start,k_end
1152	DO i = i_start,i_end
1153	! scalar was coupled with my
1154	tendency(i,k,j) = advect_tend(i,k,j) * msfty(i,j)
1155	ENDDO
1156	ENDDO
1157	ENDDO
1158
1159	DO j = j_start_spc,j_end_spc
1160	DO k = k_start_spc,k_end_spc
1161	DO i = i_start_spc,i_end_spc
1162	tendency(i,k,j) = tendency(i,k,j) + sc_tend(i,k,j,im)
1163	ENDDO
1164	ENDDO
1165	ENDDO
1166
1167	DO j = jts, min(jte,jde-1)
1168
1169	DO i = its, min(ite,ide-1)
1170	muold(i) = mu_old(i,j) + mu_base(i,j)
1171	r_munew(i) = 1./(mu_new(i,j) + mu_base(i,j))
1172	ENDDO
1173
1174	DO k = kts, min(kte,kde-1)
1175	DO i = its, min(ite,ide-1)
1176
1177	scalar_2(i,k,j,im) = (muold(i)*scalar_1(i,k,j,im) &
1178	+ dttendency(i,k,j))r_munew(i)
1179
1180	ENDDO
1181	ENDDO
1182	ENDDO
1183
1184	ENDDO
1185
1186	END IF
1187
1188	END SUBROUTINE rk_update_scalar
1189
1190	!-------------------------------------------------------------------------------
1191
1192	SUBROUTINE rk_update_scalar_pd( scs, sce, &
1193	scalar, sc_tend, &
1194	mu_old, mu_new, mu_base, &
1195	rk_step, dt, spec_zone, &
1196	config_flags, &
1197	ids, ide, jds, jde, kds, kde, &
1198	ims, ime, jms, jme, kms, kme, &
1199	its, ite, jts, jte, kts, kte )
1200
1201	IMPLICIT NONE
1202
1203	! Input data.
1204
1205	TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
1206
1207	INTEGER , INTENT(IN ) :: scs, sce, rk_step, spec_zone
1208	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1209	ims, ime, jms, jme, kms, kme, &
1210	its, ite, jts, jte, kts, kte
1211
1212	REAL, INTENT(IN ) :: dt
1213
1214	REAL, DIMENSION(ims:ime, kms:kme, jms:jme , scs:sce), &
1215	INTENT(INOUT) :: scalar, &
1216	sc_tend
1217
1218	REAL, DIMENSION(ims:ime, jms:jme ), INTENT(IN ) :: mu_old, &
1219	mu_new, &
1220	mu_base
1221
1222	INTEGER :: i,j,k,im
1223	REAL :: sc_middle, msfsq
1224	REAL, DIMENSION(its:ite) :: muold, r_munew
1225
1226	REAL, DIMENSION(its:ite, kts:kte, jts:jte ) :: tendency
1227
1228	INTEGER :: i_start,i_end,j_start,j_end,k_start,k_end
1229	INTEGER :: i_start_spc,i_end_spc,j_start_spc,j_end_spc,k_start_spc,k_end_spc
1230
1231	!<DESCRIPTION>
1232	!
1233	! rk_scalar_update advances the scalar equation given the time t value
1234	! of the scalar and the scalar tendency.
1235	!
1236	!</DESCRIPTION>
1237
1238
1239	!
1240	! set loop limits.
1241
1242	i_start = its
1243	i_end = ite
1244	j_start = jts
1245	j_end = jte
1246	k_start = kts
1247	k_end = kte-1
1248	IF(j_end == jde) j_end = j_end - 1
1249	IF(i_end == ide) i_end = i_end - 1
1250
1251	i_start_spc = i_start
1252	i_end_spc = i_end
1253	j_start_spc = j_start
1254	j_end_spc = j_end
1255	k_start_spc = k_start
1256	k_end_spc = k_end
1257
1258	IF( config_flags%nested .or. config_flags%specified ) THEN
1259	IF( .NOT. config_flags%periodic_x)i_start = max( its,ids+spec_zone )
1260	IF( .NOT. config_flags%periodic_x)i_end = min( ite,ide-spec_zone-1 )
1261	j_start = max( jts,jds+spec_zone )
1262	j_end = min( jte,jde-spec_zone-1 )
1263	k_start = kts
1264	k_end = min( kte, kde-1 )
1265	ENDIF
1266
1267	DO im = scs, sce
1268
1269	DO j = jts, min(jte,jde-1)
1270	DO k = kts, min(kte,kde-1)
1271	DO i = its, min(ite,ide-1)
1272	tendency(i,k,j) = 0.
1273	ENDDO
1274	ENDDO
1275	ENDDO
1276
1277	DO j = j_start_spc,j_end_spc
1278	DO k = k_start_spc,k_end_spc
1279	DO i = i_start_spc,i_end_spc
1280	tendency(i,k,j) = tendency(i,k,j) + sc_tend(i,k,j,im)
1281	sc_tend(i,k,j,im) = 0.
1282	ENDDO
1283	ENDDO
1284	ENDDO
1285
1286	DO j = jts, min(jte,jde-1)
1287
1288	DO i = its, min(ite,ide-1)
1289	muold(i) = mu_old(i,j) + mu_base(i,j)
1290	r_munew(i) = 1./(mu_new(i,j) + mu_base(i,j))
1291	ENDDO
1292
1293	DO k = kts, min(kte,kde-1)
1294	DO i = its, min(ite,ide-1)
1295
1296	scalar(i,k,j,im) = (muold(i)*scalar(i,k,j,im) &
1297	+ dttendency(i,k,j))r_munew(i)
1298	ENDDO
1299	ENDDO
1300	ENDDO
1301
1302	ENDDO
1303
1304	END SUBROUTINE rk_update_scalar_pd
1305
1306	!------------------------------------------------------------
1307
1308	SUBROUTINE init_zero_tendency(ru_tendf, rv_tendf, rw_tendf, ph_tendf, &
1309	t_tendf, tke_tendf, mu_tendf, &
1310	moist_tendf,chem_tendf,scalar_tendf, &
1311	n_moist,n_chem,n_scalar,rk_step, &
1312	ids, ide, jds, jde, kds, kde, &
1313	ims, ime, jms, jme, kms, kme, &
1314	its, ite, jts, jte, kts, kte )
1315	!-----------------------------------------------------------------------
1316	IMPLICIT NONE
1317	!-----------------------------------------------------------------------
1318
1319	INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1320	ims, ime, jms, jme, kms, kme, &
1321	its, ite, jts, jte, kts, kte
1322
1323	INTEGER , INTENT(IN ) :: n_moist,n_chem,n_scalar,rk_step
1324
1325	REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(INOUT) :: &
1326	ru_tendf, &
1327	rv_tendf, &
1328	rw_tendf, &
1329	ph_tendf, &
1330	t_tendf, &
1331	tke_tendf
1332
1333	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(INOUT) :: mu_tendf
1334
1335	REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_moist),INTENT(INOUT)::&
1336	moist_tendf
1337
1338	REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_chem ),INTENT(INOUT)::&
1339	chem_tendf
1340
1341	REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_scalar ),INTENT(INOUT)::&
1342	scalar_tendf
1343
1344	! LOCAL VARS
1345
1346	INTEGER :: im, ic, is
1347
1348	!<DESCRIPTION>
1349	!
1350	! init_zero_tendency
1351	! sets tendency arrays to zero for all prognostic variables.
1352	!
1353	!</DESCRIPTION>
1354
1355
1356	CALL zero_tend ( ru_tendf, &
1357	ids, ide, jds, jde, kds, kde, &
1358	ims, ime, jms, jme, kms, kme, &
1359	its, ite, jts, jte, kts, kte )
1360
1361	CALL zero_tend ( rv_tendf, &
1362	ids, ide, jds, jde, kds, kde, &
1363	ims, ime, jms, jme, kms, kme, &
1364	its, ite, jts, jte, kts, kte )
1365
1366	CALL zero_tend ( rw_tendf, &
1367	ids, ide, jds, jde, kds, kde, &
1368	ims, ime, jms, jme, kms, kme, &
1369	its, ite, jts, jte, kts, kte )
1370
1371	CALL zero_tend ( ph_tendf, &
1372	ids, ide, jds, jde, kds, kde, &
1373	ims, ime, jms, jme, kms, kme, &
1374	its, ite, jts, jte, kts, kte )
1375
1376	CALL zero_tend ( t_tendf, &
1377	ids, ide, jds, jde, kds, kde, &
1378	ims, ime, jms, jme, kms, kme, &
1379	its, ite, jts, jte, kts, kte )
1380
1381	CALL zero_tend ( tke_tendf, &
1382	ids, ide, jds, jde, kds, kde, &
1383	ims, ime, jms, jme, kms, kme, &
1384	its, ite, jts, jte, kts, kte )
1385
1386	CALL zero_tend ( mu_tendf, &
1387	ids, ide, jds, jde, kds, kds, &
1388	ims, ime, jms, jme, kms, kms, &
1389	its, ite, jts, jte, kts, kts )
1390
1391	! DO im=PARAM_FIRST_SCALAR,n_moist
1392	DO im=1,n_moist ! make sure first one is zero too
1393	CALL zero_tend ( moist_tendf(ims,kms,jms,im), &
1394	ids, ide, jds, jde, kds, kde, &
1395	ims, ime, jms, jme, kms, kme, &
1396	its, ite, jts, jte, kts, kte )
1397	ENDDO
1398
1399	! DO ic=PARAM_FIRST_SCALAR,n_chem
1400	DO ic=1,n_chem ! make sure first one is zero too
1401	CALL zero_tend ( chem_tendf(ims,kms,jms,ic), &
1402	ids, ide, jds, jde, kds, kde, &
1403	ims, ime, jms, jme, kms, kme, &
1404	its, ite, jts, jte, kts, kte )
1405	ENDDO
1406
1407	! DO ic=PARAM_FIRST_SCALAR,n_scalar
1408	DO ic=1,n_scalar ! make sure first one is zero too
1409	CALL zero_tend ( scalar_tendf(ims,kms,jms,ic), &
1410	ids, ide, jds, jde, kds, kde, &
1411	ims, ime, jms, jme, kms, kme, &
1412	its, ite, jts, jte, kts, kte )
1413	ENDDO
1414
1415	END SUBROUTINE init_zero_tendency
1416
1417	!===================================================================
1418
1419
1420	SUBROUTINE dump_data( a, field, io_unit, &
1421	ims, ime, jms, jme, kms, kme, &
1422	ids, ide, jds, jde, kds, kde )
1423	implicit none
1424	integer :: ims, ime, jms, jme, kms, kme, &
1425	ids, ide, jds, jde, kds, kde
1426	real, dimension(ims:ime, kms:kme, jds:jde) :: a
1427	character :: field
1428	integer :: io_unit
1429
1430	integer :: is,ie,js,je,ks,ke
1431
1432	!<DESCRIPTION
1433	!
1434	! quick and dirty debug io utility
1435	!
1436	!</DESCRIPTION
1437
1438	is = ids
1439	ie = ide-1
1440	js = jds
1441	je = jde-1
1442	ks = kds
1443	ke = kde-1
1444
1445	if(field == 'u') ie = ide
1446	if(field == 'v') je = jde
1447	if(field == 'w') ke = kde
1448
1449	write(io_unit) is,ie,ks,ke,js,je
1450	write(io_unit) a(is:ie, ks:ke, js:je)
1451
1452	end subroutine dump_data
1453
1454	!-----------------------------------------------------------------------
1455
1456	SUBROUTINE calculate_phy_tend (config_flags,mu,muu,muv,pi3d, &
1457	RTHRATEN, &
1458	RUBLTEN,RVBLTEN,RTHBLTEN, &
1459	RQVBLTEN,RQCBLTEN,RQIBLTEN, &
1460	RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN, &
1461	RQICUTEN,RQSCUTEN, &
1462	RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RQVNDGDTEN, &
1463	RMUNDGDTEN, &
1464	ids,ide, jds,jde, kds,kde, &
1465	ims,ime, jms,jme, kms,kme, &
1466	its,ite, jts,jte, kts,kte )
1467	!-----------------------------------------------------------------------
1468	IMPLICIT NONE
1469
1470	TYPE(grid_config_rec_type), INTENT(IN) :: config_flags
1471
1472	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1473	ims,ime, jms,jme, kms,kme, &
1474	its,ite, jts,jte, kts,kte
1475
1476	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
1477	INTENT(IN ) :: pi3d
1478
1479	REAL, DIMENSION( ims:ime, jms:jme ) , &
1480	INTENT(IN ) :: mu, &
1481	muu, &
1482	muv
1483
1484
1485	! radiation
1486
1487	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
1488	INTENT(INOUT) :: RTHRATEN
1489
1490	! cumulus
1491
1492	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
1493	INTENT(INOUT) :: &
1494	RTHCUTEN, &
1495	RQVCUTEN, &
1496	RQCCUTEN, &
1497	RQRCUTEN, &
1498	RQICUTEN, &
1499	RQSCUTEN
1500	! pbl
1501
1502	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
1503	INTENT(INOUT) :: RUBLTEN, &
1504	RVBLTEN, &
1505	RTHBLTEN, &
1506	RQVBLTEN, &
1507	RQCBLTEN, &
1508	RQIBLTEN
1509
1510	! fdda
1511
1512	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
1513	INTENT(INOUT) :: RUNDGDTEN, &
1514	RVNDGDTEN, &
1515	RTHNDGDTEN, &
1516	RQVNDGDTEN
1517	REAL, DIMENSION( ims:ime, jms:jme ) , &
1518	INTENT(INOUT) :: RMUNDGDTEN
1519
1520	INTEGER :: i,k,j
1521	INTEGER :: itf,ktf,jtf,itsu,jtsv
1522
1523	!-----------------------------------------------------------------------
1524
1525	!<DESCRIPTION>
1526	!
1527	! calculate_phy_tend couples the physics tendencies to the column mass (mu),
1528	! because prognostic equations are in flux form, but physics tendencies are
1529	! computed for uncoupled variables.
1530	!
1531	!</DESCRIPTION>
1532
1533	itf=MIN(ite,ide-1)
1534	jtf=MIN(jte,jde-1)
1535	ktf=MIN(kte,kde-1)
1536	itsu=MAX(its,ids+1)
1537	jtsv=MAX(jts,jds+1)
1538
1539	! radiation
1540
1541	IF (config_flags%ra_lw_physics .gt. 0 .or. config_flags%ra_sw_physics .gt. 0) THEN
1542
1543	DO J=jts,jtf
1544	DO K=kts,ktf
1545	DO I=its,itf
1546	RTHRATEN(I,K,J)=mu(I,J)*RTHRATEN(I,K,J)
1547	ENDDO
1548	ENDDO
1549	ENDDO
1550
1551	ENDIF
1552
1553	! cumulus
1554
1555	IF (config_flags%cu_physics .gt. 0) THEN
1556
1557	DO J=jts,jtf
1558	DO I=its,itf
1559	DO K=kts,ktf
1560	RTHCUTEN(I,K,J)=mu(I,J)*RTHCUTEN(I,K,J)
1561	RQVCUTEN(I,K,J)=mu(I,J)*RQVCUTEN(I,K,J)
1562	ENDDO
1563	ENDDO
1564	ENDDO
1565
1566	IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
1567	DO J=jts,jtf
1568	DO I=its,itf
1569	DO K=kts,ktf
1570	RQCCUTEN(I,K,J)=mu(I,J)*RQCCUTEN(I,K,J)
1571	ENDDO
1572	ENDDO
1573	ENDDO
1574	ENDIF
1575
1576	IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
1577	DO J=jts,jtf
1578	DO I=its,itf
1579	DO K=kts,ktf
1580	RQRCUTEN(I,K,J)=mu(I,J)*RQRCUTEN(I,K,J)
1581	ENDDO
1582	ENDDO
1583	ENDDO
1584	ENDIF
1585
1586	IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
1587	DO J=jts,jtf
1588	DO I=its,itf
1589	DO K=kts,ktf
1590	RQICUTEN(I,K,J)=mu(I,J)*RQICUTEN(I,K,J)
1591	ENDDO
1592	ENDDO
1593	ENDDO
1594	ENDIF
1595
1596	IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
1597	DO J=jts,jtf
1598	DO I=its,itf
1599	DO K=kts,ktf
1600	RQSCUTEN(I,K,J)=mu(I,J)*RQSCUTEN(I,K,J)
1601	ENDDO
1602	ENDDO
1603	ENDDO
1604	ENDIF
1605
1606	ENDIF
1607
1608	! pbl
1609
1610	! IF (config_flags%bl_pbl_physics .gt. 0) THEN
1611	!!****MARS
1612	IF ( (config_flags%bl_pbl_physics .gt. 0) .OR. (config_flags%modif_wrf) ) THEN
1613
1614	DO J=jts,jtf
1615	DO K=kts,ktf
1616	DO I=its,itf
1617	RUBLTEN(I,K,J) =mu(I,J)*RUBLTEN(I,K,J)
1618	RVBLTEN(I,K,J) =mu(I,J)*RVBLTEN(I,K,J)
1619	RTHBLTEN(I,K,J)=mu(I,J)*RTHBLTEN(I,K,J)
1620	ENDDO
1621	ENDDO
1622	ENDDO
1623
1624	IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
1625	DO J=jts,jtf
1626	DO K=kts,ktf
1627	DO I=its,itf
1628	RQVBLTEN(I,K,J)=mu(I,J)*RQVBLTEN(I,K,J)
1629	ENDDO
1630	ENDDO
1631	ENDDO
1632	ENDIF
1633
1634	IF (P_QC .ge. PARAM_FIRST_SCALAR) THEN
1635	DO J=jts,jtf
1636	DO K=kts,ktf
1637	DO I=its,itf
1638	RQCBLTEN(I,K,J)=mu(I,J)*RQCBLTEN(I,K,J)
1639	ENDDO
1640	ENDDO
1641	ENDDO
1642	ENDIF
1643
1644	IF (P_QI .ge. PARAM_FIRST_SCALAR) THEN
1645	DO J=jts,jtf
1646	DO K=kts,ktf
1647	DO I=its,itf
1648	RQIBLTEN(I,K,J)=mu(I,J)*RQIBLTEN(I,K,J)
1649	ENDDO
1650	ENDDO
1651	ENDDO
1652	ENDIF
1653
1654	ENDIF
1655
1656	! fdda
1657	! note fdda u and v tendencies are staggered, also only interior points have muu/muv,
1658	! so only couple those
1659
1660	IF (config_flags%grid_fdda .gt. 0) THEN
1661
1662	DO J=jts,jtf
1663	DO K=kts,ktf
1664	DO I=itsu,itf
1665	! if( i == itf/2 .AND. j == jtf/2 .AND. k == ktf/2 ) &
1666	! write(*,'(a,3i6,e15.5)') 'u_ten before=',i,k,j, RUNDGDTEN(i,k,j)
1667	RUNDGDTEN(I,K,J) =muu(I,J)*RUNDGDTEN(I,K,J)
1668	! if( i == itf/2 .AND. j == jtf/2 .AND. k==ktf/2 ) &
1669	! write(*,'(a,2f15.5)') 'mu, muu=',mu(i,j), muu(i,j)
1670	! if( i == itf/2 .AND. j == jtf/2 .AND. k == ktf/2 ) &
1671	! write(*,'(a,3i6,e15.5)') 'u_ten after=',i,k,j, RUNDGDTEN(i,k,j)
1672	! if( RUNDGDTEN(i,k,j) > 30.0 ) write(,) 'IKJ=',i,k,j
1673	ENDDO
1674	ENDDO
1675	ENDDO
1676	! write(*,'(a,e15.5)') 'u_ten MAXIMUM after=', maxval(RUNDGDTEN)
1677	DO J=jtsv,jtf
1678	DO K=kts,ktf
1679	DO I=its,itf
1680	RVNDGDTEN(I,K,J) =muv(I,J)*RVNDGDTEN(I,K,J)
1681	ENDDO
1682	ENDDO
1683	ENDDO
1684	DO J=jts,jtf
1685	DO K=kts,ktf
1686	DO I=its,itf
1687	! if( i == itf/2 .AND. j == jtf/2 .AND. k == ktf/2 ) &
1688	! write(*,'(a,3i6,e15.5)') 'th before=',i,k,j, RTHNDGDTEN(I,K,J)
1689	RTHNDGDTEN(I,K,J)=mu(I,J)*RTHNDGDTEN(I,K,J)
1690	! RMUNDGDTEN(I,J) - no coupling
1691	! if( i == itf/2 .AND. j == jtf/2 .AND. k == ktf/2 ) &
1692	! write(*,'(a,3i6,e15.5)') 'th after=',i,k,j, RTHNDGDTEN(I,K,J)
1693	ENDDO
1694	ENDDO
1695	ENDDO
1696
1697	IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
1698	DO J=jts,jtf
1699	DO K=kts,ktf
1700	DO I=its,itf
1701	RQVNDGDTEN(I,K,J)=mu(I,J)*RQVNDGDTEN(I,K,J)
1702	ENDDO
1703	ENDDO
1704	ENDDO
1705	ENDIF
1706
1707	ENDIF
1708
1709	END SUBROUTINE calculate_phy_tend
1710
1711	!-----------------------------------------------------------------------
1712
1713	SUBROUTINE positive_definite_filter ( a, &
1714	ids,ide, jds,jde, kds,kde, &
1715	ims,ime, jms,jme, kms,kme, &
1716	its,ite, jts,jte, kts,kte )
1717
1718	IMPLICIT NONE
1719
1720	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1721	ims,ime, jms,jme, kms,kme, &
1722	its,ite, jts,jte, kts,kte
1723
1724	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: a
1725
1726	INTEGER :: i,k,j
1727
1728	!<DESCRIPTION>
1729	!
1730	! debug and testing code for bounding a variable
1731	!
1732	!</DESCRIPTION>
1733
1734	DO j=jts,min(jte,jde-1)
1735	DO k=kts,kte-1
1736	DO i=its,min(ite,ide-1)
1737	! a(i,k,j) = max(a(i,k,j),0.)
1738	a(i,k,j) = min(1000.,max(a(i,k,j),0.))
1739	ENDDO
1740	ENDDO
1741	ENDDO
1742
1743	END SUBROUTINE positive_definite_filter
1744
1745	!-----------------------------------------------------------------------
1746
1747	SUBROUTINE bound_tke ( tke, tke_upper_bound, &
1748	ids,ide, jds,jde, kds,kde, &
1749	ims,ime, jms,jme, kms,kme, &
1750	its,ite, jts,jte, kts,kte )
1751
1752	IMPLICIT NONE
1753
1754	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1755	ims,ime, jms,jme, kms,kme, &
1756	its,ite, jts,jte, kts,kte
1757
1758	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: tke
1759	REAL, INTENT( IN) :: tke_upper_bound
1760
1761	INTEGER :: i,k,j
1762
1763	!<DESCRIPTION>
1764	!
1765	! bounds tke between zero and tke_upper_bound.
1766	!
1767	!</DESCRIPTION>
1768
1769	DO j=jts,min(jte,jde-1)
1770	DO k=kts,kte-1
1771	DO i=its,min(ite,ide-1)
1772	tke(i,k,j) = min(tke_upper_bound,max(tke(i,k,j),0.))
1773	ENDDO
1774	ENDDO
1775	ENDDO
1776
1777	END SUBROUTINE bound_tke
1778
1779
1780
1781	END MODULE module_em

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

Download in other formats:

Original Format