Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_em.F @ 3094

Last change on this file since 3094 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

Rev	Line
[17]	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