Context Navigation

solve_em.F @ 200

Last change on this file since 200 was 184, checked in by lfita, 11 years ago

Adding 'qc' on mp==0 in Registry/Registry?.EM.LMDZ. With this now we have also

'QCLOUD' even with mp=0. This is the way to keep liquid water for the next
time-step and have large scale precipitation!

Adding some checks around moisture values and tendencies...

File size: 226.1 KB

Line
1	!WRF:MEDIATION_LAYER:SOLVER
2
3	SUBROUTINE solve_em ( grid , config_flags &
4	! Arguments generated from Registry
5	#include "dummy_new_args.inc"
6	!
7	)
8	! Driver layer modules
9	USE module_state_description
10	USE module_domain, ONLY : &
11	domain, get_ijk_from_grid, get_ijk_from_subgrid &
12	,domain_get_current_time, domain_get_start_time &
13	,domain_get_sim_start_time, domain_clock_get
14	USE module_domain_type, ONLY : history_alarm
15	USE module_configure, ONLY : grid_config_rec_type
16	USE module_driver_constants
17	USE module_machine
18	USE module_tiles, ONLY : set_tiles
19	#ifdef DM_PARALLEL
20	USE module_dm, ONLY : &
21	local_communicator, mytask, ntasks, ntasks_x, ntasks_y &
22	,local_communicator_periodic, wrf_dm_maxval
23	USE module_comm_dm, ONLY : &
24	halo_em_a_sub,halo_em_b_sub,halo_em_c2_sub,halo_em_chem_e_3_sub &
25	,halo_em_chem_e_5_sub,halo_em_chem_e_7_sub,halo_em_chem_old_e_5_sub &
26	,halo_em_chem_old_e_7_sub,halo_em_c_sub,halo_em_d2_3_sub &
27	,halo_em_d2_5_sub,halo_em_d3_3_sub,halo_em_d3_5_sub,halo_em_d_sub &
28	,halo_em_e_3_sub,halo_em_e_5_sub,halo_em_hydro_uv_sub &
29	,halo_em_moist_e_3_sub,halo_em_moist_e_5_sub,halo_em_moist_e_7_sub &
30	,halo_em_moist_old_e_5_sub,halo_em_moist_old_e_7_sub &
31	,halo_em_scalar_e_3_sub,halo_em_scalar_e_5_sub,halo_em_scalar_e_7_sub &
32	,halo_em_scalar_old_e_5_sub,halo_em_scalar_old_e_7_sub,halo_em_tke_3_sub &
33	,halo_em_tke_5_sub,halo_em_tke_7_sub,halo_em_tke_advect_3_sub &
34	,halo_em_tke_advect_5_sub,halo_em_tke_old_e_5_sub &
35	,halo_em_tke_old_e_7_sub,halo_em_tracer_e_3_sub,halo_em_tracer_e_5_sub &
36	,halo_em_tracer_e_7_sub,halo_em_tracer_old_e_5_sub &
37	,halo_em_tracer_old_e_7_sub,period_bdy_em_a_sub &
38	,period_bdy_em_b3_sub,period_bdy_em_b_sub,period_bdy_em_chem2_sub &
39	,period_bdy_em_chem_old_sub,period_bdy_em_chem_sub,period_bdy_em_d3_sub &
40	,period_bdy_em_d_sub,period_bdy_em_e_sub,period_bdy_em_moist2_sub &
41	,period_bdy_em_moist_old_sub,period_bdy_em_moist_sub &
42	,period_bdy_em_scalar2_sub,period_bdy_em_scalar_old_sub &
43	,period_bdy_em_scalar_sub,period_bdy_em_tke_old_sub &
44	,period_bdy_em_tracer2_sub,period_bdy_em_tracer_old_sub &
45	,period_bdy_em_tracer_sub,period_em_da_sub,period_em_hydro_uv_sub
46	#endif
47	USE module_utility
48	! Mediation layer modules
49	! Model layer modules
50	USE module_model_constants
51	USE module_small_step_em
52	USE module_em
53	USE module_big_step_utilities_em
54	USE module_bc
55	USE module_bc_em
56	USE module_solvedebug_em
57	USE module_physics_addtendc
58	USE module_diffusion_em
59	USE module_polarfft
60	USE module_microphysics_driver
61	USE module_microphysics_zero_out
62	USE module_fddaobs_driver
63	USE module_diagnostics
64	#ifdef WRF_CHEM
65	USE module_input_chem_data
66	USE module_input_tracer
67	USE module_chem_utilities
68	#endif
69	USE module_first_rk_step_part1
70	USE module_first_rk_step_part2
71	USE module_llxy, ONLY : proj_cassini
72	USE module_avgflx_em, ONLY : zero_avgflx, upd_avgflx
73	#ifdef LMDZ
74	USE module_lmdz_phys
75	USE module_domain_type, ONLY: RESTART_ALARM
76	USE module_domain, ONLY: domain_get_time_since_sim_start
77	USE module_streams
78
79	#endif
80
81	IMPLICIT NONE
82
83	! Input data.
84
85	TYPE(domain) , TARGET :: grid
86
87	! Definitions of dummy arguments to this routine (generated from Registry).
88	#include "dummy_new_decl.inc"
89
90	! Structure that contains run-time configuration (namelist) data for domain
91	TYPE (grid_config_rec_type) , INTENT(IN) :: config_flags
92
93	! Local data
94
95	INTEGER :: k_start , k_end, its, ite, jts, jte
96	INTEGER :: ids , ide , jds , jde , kds , kde , &
97	ims , ime , jms , jme , kms , kme , &
98	ips , ipe , jps , jpe , kps , kpe
99
100	INTEGER :: sids , side , sjds , sjde , skds , skde , &
101	sims , sime , sjms , sjme , skms , skme , &
102	sips , sipe , sjps , sjpe , skps , skpe
103
104
105	INTEGER :: imsx, imex, jmsx, jmex, kmsx, kmex, &
106	ipsx, ipex, jpsx, jpex, kpsx, kpex, &
107	imsy, imey, jmsy, jmey, kmsy, kmey, &
108	ipsy, ipey, jpsy, jpey, kpsy, kpey
109
110	INTEGER :: ij , iteration
111	INTEGER :: im , num_3d_m , ic , num_3d_c , is , num_3d_s
112	INTEGER :: loop
113	INTEGER :: sz
114	INTEGER :: iswater
115
116	LOGICAL :: specified_bdy, channel_bdy
117
118	REAL :: t_new
119
120	! Changes in tendency at this timestep
121	real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,grid%sm33:grid%em33) :: h_tendency, &
122	z_tendency
123
124	! Whether advection should produce decoupled horizontal and vertical advective tendency outputs
125	LOGICAL :: tenddec
126
127	#ifdef WRF_CHEM
128	! Index cross-referencing array for tendency accumulation
129	INTEGER, DIMENSION( num_chem ) :: adv_ct_indices
130	#endif
131
132	! storage for tendencies and decoupled state (generated from Registry)
133
134	#include <i1_decl.inc>
135	! Previous time level of tracer arrays now defined as i1 variables;
136	! the state 4d arrays now redefined as 1-time level arrays in Registry.
137	! Benefit: save memory in nested runs, since only 1 domain is active at a
138	! time. Potential problem on stack-limited architectures: increases
139	! amount of data on program stack by making these automatic arrays.
140
141	INTEGER :: rc
142	INTEGER :: number_of_small_timesteps, rk_step
143	INTEGER :: klevel,ijm,ijp,i,j,k,size1,size2 ! for prints/plots only
144	INTEGER :: idum1, idum2, dynamics_option
145
146	INTEGER :: rk_order, iwmax, jwmax, kwmax
147	REAL :: dt_rk, dts_rk, dts, dtm, wmax
148	REAL , ALLOCATABLE , DIMENSION(:) :: max_vert_cfl_tmp, max_horiz_cfl_tmp
149	LOGICAL :: leapfrog
150	INTEGER :: l,kte,kk
151	LOGICAL :: f_flux ! flag for computing averaged fluxes in cu_gd
152	REAL :: curr_secs
153	INTEGER :: num_sound_steps
154	INTEGER :: idex, jdex
155	REAL :: max_msft
156	REAL :: spacing
157
158	INTEGER :: ii, jj !kk is above after l,kte
159	REAL :: dclat
160	INTEGER :: debug_level
161
162	! urban related variables
163	INTEGER :: NUM_ROOF_LAYERS, NUM_WALL_LAYERS, NUM_ROAD_LAYERS ! urban
164
165	TYPE(WRFU_TimeInterval) :: tmpTimeInterval
166	REAL :: real_time
167	LOGICAL :: adapt_step_flag
168	LOGICAL :: fill_w_flag
169
170	! variables for flux-averaging code 20091223
171	CHARACTER*256 :: message, message2
172	REAL :: old_dt
173	TYPE(WRFU_Time) :: temp_time, CurrTime
174	INTEGER, PARAMETER :: precision = 100
175	INTEGER :: num, den
176	TYPE(WRFU_TimeInterval) :: dtInterval
177
178	#ifdef LMDZ
179	INTEGER :: im2, jm2, km2
180	INTEGER :: ix,iy,iz
181	INTEGER :: days, seconds, Sn, Sd
182	TYPE(WRFU_TimeInterval) :: timeSinceSimStart
183	TYPE(WRFU_Time) :: initime, simtime
184	LOGICAL :: wrftestrst, wrftestin
185	REAL :: minSinceSimStart
186	CHARACTER(LEN=256) :: mminlu, mminsl
187	CHARACTER(LEN=50) :: errmsg
188	INTEGER :: hr, minute, sec, ms, julyr, &
189	julday
190	REAL :: gmt
191
192
193	errmsg = 'ERROR -- error -- ERROR -- error'
194
195	im2 = config_flags%i_check_point
196	jm2 = config_flags%j_check_point
197	km2 = config_flags%k_check_point
198	#endif
199
200	! Define benchmarking timers if -DBENCH is compiled
201	#include <bench_solve_em_def.h>
202
203	!----------------------
204	! Executable statements
205	!----------------------
206
207	!<DESCRIPTION>
208	!<pre>
209	! solve_em is the main driver for advancing a grid a single timestep.
210	! It is a mediation-layer routine -> DM and SM calls are made where
211	! needed for parallel processing.
212	!
213	! solve_em can integrate the equations using 3 time-integration methods
214	!
215	! - 3rd order Runge-Kutta time integration (recommended)
216	!
217	! - 2nd order Runge-Kutta time integration
218	!
219	! The main sections of solve_em are
220	!
221	! (1) Runge-Kutta (RK) loop
222	!
223	! (2) Non-timesplit physics (i.e., tendencies computed for updating
224	! model state variables during the first RK sub-step (loop)
225	!
226	! (3) Small (acoustic, sound) timestep loop - within the RK sub-steps
227	!
228	! (4) scalar advance for moist and chem scalar variables (and TKE)
229	! within the RK sub-steps.
230	!
231	! (5) time-split physics (after the RK step), currently this includes
232	! only microphyics
233	!
234	! A more detailed description of these sections follows.
235	!</pre>
236	!</DESCRIPTION>
237
238	! Initialize timers if compiled with -DBENCH
239	#include <bench_solve_em_init.h>
240
241	! set runge-kutta solver (2nd or 3rd order)
242
243	dynamics_option = config_flags%rk_ord
244
245	! Obtain dimension information stored in the grid data structure.
246
247	CALL get_ijk_from_grid ( grid , &
248	ids, ide, jds, jde, kds, kde, &
249	ims, ime, jms, jme, kms, kme, &
250	ips, ipe, jps, jpe, kps, kpe, &
251	imsx, imex, jmsx, jmex, kmsx, kmex, &
252	ipsx, ipex, jpsx, jpex, kpsx, kpex, &
253	imsy, imey, jmsy, jmey, kmsy, kmey, &
254	ipsy, ipey, jpsy, jpey, kpsy, kpey )
255
256	CALL get_ijk_from_subgrid ( grid , &
257	sids, side, sjds, sjde, skds, skde, &
258	sims, sime, sjms, sjme, skms, skme, &
259	sips, sipe, sjps, sjpe, skps, skpe )
260	k_start = kps
261	k_end = kpe
262
263	num_3d_m = num_moist
264	num_3d_c = num_chem
265	num_3d_s = num_scalar
266
267	f_flux = config_flags%do_avgflx_cugd .EQ. 1
268
269	! Compute these starting and stopping locations for each tile and number of tiles.
270	! See: http://www.mmm.ucar.edu/wrf/WG2/topics/settiles
271	CALL set_tiles ( grid , ids , ide , jds , jde , ips , ipe , jps , jpe )
272
273	! Max values of CFL for adaptive time step scheme
274
275	ALLOCATE (max_vert_cfl_tmp(grid%num_tiles))
276	ALLOCATE (max_horiz_cfl_tmp(grid%num_tiles))
277
278	!
279	! Calculate current time in seconds since beginning of model run.
280	! Unfortunately, ESMF does not seem to have a way to return
281	! floating point seconds based on a TimeInterval. So, we will
282	! calculate it here--but, this is not clean!!
283	!
284	tmpTimeInterval = domain_get_current_time ( grid ) - domain_get_sim_start_time ( grid )
285	curr_secs = real_time(tmpTimeInterval)
286
287	old_dt = grid%dt ! store old time step for flux averaging code at end of RK loop
288	!-----------------------------------------------------------------------------
289	! Adaptive time step: Added by T. Hutchinson, WSI 3/5/07
290	! In this call, we do the time-step adaptation and set time-dependent lateral
291	! boundary condition nudging weights.
292	!
293	IF ( (config_flags%use_adaptive_time_step) .and. &
294	( (.not. grid%nested) .or. &
295	( (grid%nested) .and. (abs(grid%dtbc) < 0.0001) ) ) )THEN
296	CALL adapt_timestep(grid, config_flags)
297	adapt_step_flag = .TRUE.
298	ELSE
299	adapt_step_flag = .FALSE.
300	ENDIF
301	! End of adaptive time step modifications
302	!-----------------------------------------------------------------------------
303
304	grid%itimestep = grid%itimestep + 1
305
306	IF (config_flags%polar) dclat = 90./REAL(jde-jds) !(0.5 * 180/ny)
307
308	#ifdef WRF_CHEM
309
310	kte=min(k_end,kde-1)
311	# ifdef DM_PARALLEL
312	if ( num_chem >= PARAM_FIRST_SCALAR ) then
313	!-----------------------------------------------------------------------
314	! see matching halo calls below for stencils
315	!--------------------------------------------------------------
316	CALL wrf_debug ( 200 , ' call HALO_RK_CHEM' )
317	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
318	# include "HALO_EM_CHEM_E_3.inc"
319	IF( config_flags%progn > 0 ) THEN
320	# include "HALO_EM_SCALAR_E_3.inc"
321	ENDIF
322	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
323	# include "HALO_EM_CHEM_E_5.inc"
324	IF( config_flags%progn > 0 ) THEN
325	# include "HALO_EM_SCALAR_E_5.inc"
326	ENDIF
327	ELSE
328	WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order
329	CALL wrf_error_fatal(TRIM(wrf_err_message))
330	ENDIF
331	ENDIF
332	if ( num_tracer >= PARAM_FIRST_SCALAR ) then
333	!-----------------------------------------------------------------------
334	! see matching halo calls below for stencils
335	!--------------------------------------------------------------
336	CALL wrf_debug ( 200 , ' call HALO_RK_tracer' )
337	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
338	# include "HALO_EM_TRACER_E_3.inc"
339	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
340	# include "HALO_EM_TRACER_E_5.inc"
341	ELSE
342	WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order
343	CALL wrf_error_fatal(TRIM(wrf_err_message))
344	ENDIF
345	ENDIF
346	# endif
347	!--------------------------------------------------------------
348	adv_ct_indices( : ) = 1
349	IF ( config_flags%chemdiag == USECHEMDIAG ) THEN
350	! modify tendency list here
351	! note that the referencing direction here is opposite of that in chem_driver
352	adv_ct_indices(p_co ) = p_advh_co
353	adv_ct_indices(p_o3 ) = p_advh_o3
354	adv_ct_indices(p_no ) = p_advh_no
355	adv_ct_indices(p_no2 ) = p_advh_no2
356	adv_ct_indices(p_hno3) = p_advh_hno3
357	adv_ct_indices(p_iso ) = p_advh_iso
358	adv_ct_indices(p_ho ) = p_advh_ho
359	adv_ct_indices(p_ho2 ) = p_advh_ho2
360	END IF
361	#endif
362
363	rk_order = config_flags%rk_ord
364
365	IF ( grid%time_step_sound == 0 ) THEN
366	! This function will give 4 for 6dx and 6 for 10dx and returns even numbers only
367	spacing = min(grid%dx, grid%dy)
368	IF ( ( config_flags%use_adaptive_time_step ) .AND. ( config_flags%map_proj == PROJ_CASSINI ) ) THEN
369	max_msft=MIN ( MAX(grid%max_msftx, grid%max_msfty) , &
370	1.0/COS(config_flags%fft_filter_lat*degrad) )
371	num_sound_steps = max ( 2 * ( INT (300. * grid%dt / (spacing / max_msft) - 0.01 ) + 1 ), 4 )
372	ELSE IF ( config_flags%use_adaptive_time_step ) THEN
373	max_msft= MAX(grid%max_msftx, grid%max_msfty)
374	num_sound_steps = max ( 2 * ( INT (300. * grid%dt / (spacing / max_msft) - 0.01 ) + 1 ), 4 )
375	ELSE
376	num_sound_steps = max ( 2 * ( INT (300. * grid%dt / spacing - 0.01 ) + 1 ), 4 )
377	END IF
378	WRITE(wrf_err_message,*)'grid spacing, dt, time_step_sound=',spacing,grid%dt,num_sound_steps
379	CALL wrf_debug ( 50 , wrf_err_message )
380	ELSE
381	num_sound_steps = grid%time_step_sound
382	ENDIF
383
384	dts = grid%dt/float(num_sound_steps)
385
386	IF (config_flags%use_adaptive_time_step) THEN
387
388	CALL get_wrf_debug_level( debug_level )
389	IF ((config_flags%time_step < 0) .AND. (debug_level.GE.50)) THEN
390	#ifdef DM_PARALLEL
391	CALL wrf_dm_maxval(grid%max_vert_cfl, idex, jdex)
392	#endif
393	WRITE(wrf_err_message,*)'variable dt, max horiz cfl, max vert cfl: ',&
394	grid%dt, grid%max_horiz_cfl, grid%max_vert_cfl
395	CALL wrf_debug ( 0 , wrf_err_message )
396	ENDIF
397
398	grid%max_cfl_val = 0
399	grid%max_horiz_cfl = 0
400	grid%max_vert_cfl = 0
401	ENDIF
402
403	! setting bdy tendencies to zero for DFI if constant_bc = true
404
405	!$OMP PARALLEL DO &
406	!$OMP PRIVATE ( ij )
407	DO ij = 1 , grid%num_tiles
408
409	! IF( config_flags%specified .AND. grid%dfi_opt .NE. DFI_NODFI &
410	! .AND. config_flags%constant_bc .AND. (grid%dfi_stage .EQ. DFI_BCK .OR. grid%dfi_stage .EQ. DFI_FWD) ) THEN
411	IF( config_flags%specified .AND. config_flags%constant_bc ) THEN
412
413	CALL zero_bdytend (grid%u_btxs,grid%u_btxe,grid%u_btys,grid%u_btye, &
414	grid%v_btxs,grid%v_btxe,grid%v_btys,grid%v_btye, &
415	grid%ph_btxs,grid%ph_btxe,grid%ph_btys,grid%ph_btye, &
416	grid%t_btxs,grid%t_btxe,grid%t_btys,grid%t_btye, &
417	grid%w_btxs,grid%w_btxe,grid%w_btys,grid%w_btye, &
418	grid%mu_btxs,grid%mu_btxe,grid%mu_btys,grid%mu_btye, &
419	moist_btxs,moist_btxe, &
420	moist_btys,moist_btye, &
421	grid%spec_bdy_width,num_3d_m, &
422	ids,ide, jds,jde, kds,kde, &
423	ims,ime, jms,jme, kms,kme, &
424	ips,ipe, jps,jpe, kps,kpe, &
425	grid%i_start(ij), grid%i_end(ij), &
426	grid%j_start(ij), grid%j_end(ij), &
427	k_start, k_end )
428
429	ENDIF
430	ENDDO
431	!$OMP END PARALLEL DO
432
433	!**********************************************************************
434	!
435	! LET US BEGIN.......
436	!
437	!<DESCRIPTION>
438	!<pre>
439	! (1) RK integration loop is named the "Runge_Kutta_loop:"
440	!
441	! Predictor-corrector type time integration.
442	! Advection terms are evaluated at time t for the predictor step,
443	! and advection is re-evaluated with the latest predicted value for
444	! each succeeding time corrector step
445	!
446	! 2nd order Runge Kutta (rk_order = 2):
447	! Step 1 is taken to the midpoint predictor, step 2 is the full step.
448	!
449	! 3rd order Runge Kutta (rk_order = 3):
450	! Step 1 is taken to from t to dt/3, step 2 is from t to dt/2,
451	! and step 3 is from t to dt.
452	!
453	! non-timesplit physics are evaluated during first RK step and
454	! these physics tendencies are stored for use in each RK pass.
455	!</pre>
456	!</DESCRIPTION>
457	!**********************************************************************
458
459	Runge_Kutta_loop: DO rk_step = 1, rk_order
460
461	! Set the step size and number of small timesteps for
462	! each part of the timestep
463
464	dtm = grid%dt
465	IF ( rk_order == 1 ) THEN
466
467	write(wrf_err_message,*)' leapfrog removed, error exit for dynamics_option = ',dynamics_option
468	CALL wrf_error_fatal( wrf_err_message )
469
470	ELSE IF ( rk_order == 2 ) THEN ! 2nd order Runge-Kutta timestep
471
472	IF ( rk_step == 1) THEN
473	dt_rk = 0.5*grid%dt
474	dts_rk = dts
475	number_of_small_timesteps = num_sound_steps/2
476	ELSE
477	dt_rk = grid%dt
478	dts_rk = dts
479	number_of_small_timesteps = num_sound_steps
480	ENDIF
481
482	ELSE IF ( rk_order == 3 ) THEN ! third order Runge-Kutta
483
484	IF ( rk_step == 1) THEN
485	dt_rk = grid%dt/3.
486	dts_rk = dt_rk
487	number_of_small_timesteps = 1
488	ELSE IF (rk_step == 2) THEN
489	dt_rk = 0.5*grid%dt
490	dts_rk = dts
491	number_of_small_timesteps = num_sound_steps/2
492	ELSE
493	dt_rk = grid%dt
494	dts_rk = dts
495	number_of_small_timesteps = num_sound_steps
496	ENDIF
497
498	ELSE
499
500	write(wrf_err_message,*)' unknown solver, error exit for dynamics_option = ',dynamics_option
501	CALL wrf_error_fatal( wrf_err_message )
502
503	END IF
504
505	! Ensure that polar meridional velocity is zero
506	IF (config_flags%polar) THEN
507	!$OMP PARALLEL DO &
508	!$OMP PRIVATE ( ij )
509	DO ij = 1 , grid%num_tiles
510	CALL zero_pole ( grid%v_1, &
511	ids, ide, jds, jde, kds, kde, &
512	ims, ime, jms, jme, kms, kme, &
513	grid%i_start(ij), grid%i_end(ij), &
514	grid%j_start(ij), grid%j_end(ij), &
515	k_start, k_end )
516	CALL zero_pole ( grid%v_2, &
517	ids, ide, jds, jde, kds, kde, &
518	ims, ime, jms, jme, kms, kme, &
519	grid%i_start(ij), grid%i_end(ij), &
520	grid%j_start(ij), grid%j_end(ij), &
521	k_start, k_end )
522	END DO
523	!$OMP END PARALLEL DO
524	END IF
525	!
526	! Time level t is in the *_2 variable in the first part
527	! of the step, and in the *_1 variable after the predictor.
528	! the latest predicted values are stored in the *_2 variables.
529	!
530	#ifdef LMDZ1
531	WRITE(message, *)' dyn_em: before rk_step_prep rk_step=', rk_step
532	CALL wrf_debug(200, message)
533	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
534	' u_tend: ', ru_tendf(im2,1,jm2)
535	CALL wrf_debug(200, message)
536	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
537	'p sfc: ',p8w(im2,kms,jm2)
538	CALL wrf_debug(200, message)
539	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
540	CALL wrf_debug(200, message)
541	#endif
542
543	CALL wrf_debug ( 200 , ' call rk_step_prep ' )
544
545	BENCH_START(step_prep_tim)
546	!$OMP PARALLEL DO &
547	!$OMP PRIVATE ( ij )
548
549	DO ij = 1 , grid%num_tiles
550
551	CALL rk_step_prep ( config_flags, rk_step, &
552	grid%u_2, grid%v_2, grid%w_2, grid%t_2, grid%ph_2, grid%mu_2, &
553	moist, &
554	grid%ru, grid%rv, grid%rw, grid%ww, grid%php, grid%alt, grid%muu, grid%muv, &
555	grid%mub, grid%mut, grid%phb, grid%pb, grid%p, grid%al, grid%alb, &
556	cqu, cqv, cqw, &
557	grid%msfux, grid%msfuy, grid%msfvx, grid%msfvx_inv, &
558	grid%msfvy, grid%msftx, grid%msfty, &
559	grid%fnm, grid%fnp, grid%dnw, grid%rdx, grid%rdy, &
560	num_3d_m, &
561	ids, ide, jds, jde, kds, kde, &
562	ims, ime, jms, jme, kms, kme, &
563	grid%i_start(ij), grid%i_end(ij), &
564	grid%j_start(ij), grid%j_end(ij), &
565	k_start, k_end )
566
567	END DO
568	!$OMP END PARALLEL DO
569	BENCH_END(step_prep_tim)
570
571	#ifdef LMDZ1
572	WRITE(message, *)' dyn_em: after rk_step_prep'
573	CALL wrf_debug(200, message)
574	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
575	' u_tend: ', ru_tendf(im2,1,jm2)
576	CALL wrf_debug(200, message)
577	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
578	'p sfc: ',p8w(im2,kms,jm2)
579	CALL wrf_debug(200, message)
580	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
581	CALL wrf_debug(200, message)
582	#endif
583
584	#ifdef DM_PARALLEL
585	!-----------------------------------------------------------------------
586	! Stencils for patch communications (WCS, 29 June 2001)
587	! Note: the small size of this halo exchange reflects the
588	! fact that we are carrying the uncoupled variables
589	! as state variables in the mass coordinate model, as
590	! opposed to the coupled variables as in the height
591	! coordinate model.
592	!
593	! * * * * *
594	! * * * * * * * * *
595	! * + * * + * * * + * *
596	! * * * * * * * * *
597	! * * * * *
598	!
599	! 3D variables - note staggering! ru(X), rv(Y), ww(Z), php(Z)
600	!
601	! ru x
602	! rv x
603	! ww x
604	! php x
605	! alt x
606	! ph_2 x
607	! phb x
608	!
609	! the following are 2D (xy) variables
610	!
611	! muu x
612	! muv x
613	! mut x
614	!--------------------------------------------------------------
615	# include "HALO_EM_A.inc"
616	#endif
617
618	! set boundary conditions on variables
619	! from big_step_prep for use in big_step_proc
620
621	#ifdef DM_PARALLEL
622	# include "PERIOD_BDY_EM_A.inc"
623	#endif
624
625	BENCH_START(set_phys_bc_tim)
626	!$OMP PARALLEL DO &
627	!$OMP PRIVATE ( ij, ii, jj, kk )
628
629	DO ij = 1 , grid%num_tiles
630
631	CALL wrf_debug ( 200 , ' call rk_phys_bc_dry_1' )
632
633	CALL rk_phys_bc_dry_1( config_flags, grid%ru, grid%rv, grid%rw, grid%ww, &
634	grid%muu, grid%muv, grid%mut, grid%php, grid%alt, grid%p, &
635	ids, ide, jds, jde, kds, kde, &
636	ims, ime, jms, jme, kms, kme, &
637	ips, ipe, jps, jpe, kps, kpe, &
638	grid%i_start(ij), grid%i_end(ij), &
639	grid%j_start(ij), grid%j_end(ij), &
640	k_start, k_end )
641	CALL set_physical_bc3d( grid%al, 'p', config_flags, &
642	ids, ide, jds, jde, kds, kde, &
643	ims, ime, jms, jme, kms, kme, &
644	ips, ipe, jps, jpe, kps, kpe, &
645	grid%i_start(ij), grid%i_end(ij), &
646	grid%j_start(ij), grid%j_end(ij), &
647	k_start , k_end )
648	CALL set_physical_bc3d( grid%ph_2, 'w', config_flags, &
649	ids, ide, jds, jde, kds, kde, &
650	ims, ime, jms, jme, kms, kme, &
651	ips, ipe, jps, jpe, kps, kpe, &
652	grid%i_start(ij), grid%i_end(ij), &
653	grid%j_start(ij), grid%j_end(ij), &
654	k_start, k_end )
655
656	#ifdef LMDZ1
657	WRITE(message, *)' dyn_em: before polar'
658	CALL wrf_debug(200, message)
659	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
660	' u_tend: ', ru_tendf(im2,1,jm2)
661	CALL wrf_debug(200, message)
662	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
663	'p sfc: ',p8w(im2,kms,jm2)
664	CALL wrf_debug(200, message)
665	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
666	CALL wrf_debug(200, message)
667	#endif
668
669	IF (config_flags%polar) THEN
670
671	!-------------------------------------------------------
672	! lat-lon grid pole-point (v) specification (extrapolate v, rv to the pole)
673	!-------------------------------------------------------
674
675	CALL pole_point_bc ( grid%v_1, &
676	ids, ide, jds, jde, kds, kde, &
677	ims, ime, jms, jme, kms, kme, &
678	grid%i_start(ij), grid%i_end(ij), &
679	grid%j_start(ij), grid%j_end(ij), &
680	k_start, k_end )
681
682	CALL pole_point_bc ( grid%v_2, &
683	ids, ide, jds, jde, kds, kde, &
684	ims, ime, jms, jme, kms, kme, &
685	grid%i_start(ij), grid%i_end(ij), &
686	grid%j_start(ij), grid%j_end(ij), &
687	k_start, k_end )
688
689	!-------------------------------------------------------
690	! end lat-lon grid pole-point (v) specification
691	!-------------------------------------------------------
692
693	ENDIF
694	END DO
695	!$OMP END PARALLEL DO
696
697	BENCH_END(set_phys_bc_tim)
698
699	rk_step_is_one : IF (rk_step == 1) THEN ! only need to initialize diffusion tendencies
700
701	!<DESCRIPTION>
702	!<pre>
703	!(2) The non-timesplit physics begins with a call to "phy_prep"
704	! (which computes some diagnostic variables such as temperature,
705	! pressure, u and v at p points, etc). This is followed by
706	! calls to the physics drivers:
707	!
708	! radiation,
709	! surface,
710	! pbl,
711	! cumulus,
712	! fddagd,
713	! 3D TKE and mixing.
714	!<pre>
715	!</DESCRIPTION>
716
717	#ifdef LMDZ
718	WRITE(message, *)' dyn_em: pre step1'
719	CALL wrf_debug(200, message)
720	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
721	' u_tend: ', ru_tendf(im2,1,jm2)
722	CALL wrf_debug(200, message)
723	IF (config_flags%lmdz_physics) THEN
724	IF (config_flags%mp_physics + config_flags%ra_lw_physics + &
725	config_flags%ra_sw_physics + config_flags%sf_sfclay_physics + &
726	config_flags%bl_pbl_physics + config_flags%cu_physics /= 0) THEN
727	PRINT *,TRIM(errmsg)
728	PRINT *,' LMDZ physics are selected from namelist. lmdz_physics= ', &
729	config_flags%lmdz_physics
730	PRINT *,' Which requires no WRF physics schemes [0 value] and tey are not:'
731	PRINT '(2(A12,1x,I2,8x))',' mp =', config_flags%mp_physics,'ra_lw =', &
732	config_flags%ra_lw_physics
733	PRINT '(2(A12,1x,I2,8x))',' ra_sw =', config_flags%ra_sw_physics, &
734	'sf_sfclay =',config_flags%sf_sfclay_physics
735	PRINT '(2(A12,1x,I2,8x))',' bl_pbl =', config_flags%bl_pbl_physics, &
736	'cu =', config_flags%cu_physics
737	message = 'WRONG namelist ste-up'
738	CALL wrf_error_fatal(TRIM(message))
739	END IF
740
741	END IF
742	WRITE(message, *)' dyn_em: before step1'
743	CALL wrf_debug(200, message)
744	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
745	' u_tend: ', ru_tendf(im2,1,jm2)
746	CALL wrf_debug(200, message)
747	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
748	'p sfc: ',p8w(im2,kms,jm2)
749	CALL wrf_debug(200, message)
750	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
751	CALL wrf_debug(200, message)
752
753	#endif
754
755	CALL first_rk_step_part1 ( grid, config_flags &
756	, moist , moist_tend &
757	, chem , chem_tend &
758	, tracer, tracer_tend &
759	, scalar , scalar_tend &
760	, fdda3d, fdda2d &
761	, ru_tendf, rv_tendf &
762	, rw_tendf, t_tendf &
763	, ph_tendf, mu_tendf &
764	, tke_tend &
765	, adapt_step_flag , curr_secs &
766	, psim , psih , wspd , gz1oz0 &
767	, br , chklowq &
768	, cu_act_flag , hol , th_phy &
769	, pi_phy , p_phy , grid%t_phy &
770	, u_phy , v_phy &
771	, dz8w , p8w , t8w , rho_phy , rho &
772	, ids, ide, jds, jde, kds, kde &
773	, ims, ime, jms, jme, kms, kme &
774	, ips, ipe, jps, jpe, kps, kpe &
775	, imsx, imex, jmsx, jmex, kmsx, kmex &
776	, ipsx, ipex, jpsx, jpex, kpsx, kpex &
777	, imsy, imey, jmsy, jmey, kmsy, kmey &
778	, ipsy, ipey, jpsy, jpey, kpsy, kpey &
779	, k_start , k_end &
780	, f_flux &
781	)
782
783	#ifdef LMDZ1
784	WRITE(message, *)' dyn_em: post step1 pre lmdz'
785	CALL wrf_debug(200, message)
786	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
787	' u_tend: ', ru_tendf(im2,1,jm2)
788	CALL wrf_debug(200, message)
789	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
790	'p sfc: ',p8w(im2,kms,jm2)
791	CALL wrf_debug(200, message)
792	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
793	CALL wrf_debug(200, message)
794	#endif
795
796	#ifdef DM_PARALLEL
797	IF ( config_flags%bl_pbl_physics == MYNNPBLSCHEME2 .OR. &
798	config_flags%bl_pbl_physics == MYNNPBLSCHEME3 ) THEN
799	# include "HALO_EM_SCALAR_E_5.inc"
800	ENDIF
801	#endif
802
803	#ifdef LMDZ
804	!!! Using grid%clock & domain_get_sim_start_time
805	!!!!
806	initime = domain_get_start_time(grid)
807	CALL domain_clock_get(grid, current_time = simtime)
808	!! timeSinceSimStart = domain_get_time_since_sim_start( grid )
809	timeSinceSimStart = simtime - initime
810	CALL WRFU_TimeIntervalGet( timeSinceSimStart, &
811	D=days, S=seconds, Sn=Sn, Sd=Sd, rc=rc )
812	IF ( rc /= WRFU_SUCCESS ) THEN
813	CALL wrf_error_fatal ( &
814	'domain_clock_get: WRFU_TimeIntervalGet() failed' )
815	ENDIF
816	! get rid of hard-coded constants
817	minSinceSimStart = ( REAL( days ) * 24. * 60. ) + &
818	( REAL( seconds ) / 60. )
819	IF ( Sd /= 0 ) THEN
820	minSinceSimStart = minSinceSimStart + &
821	( ( REAL( Sn ) / REAL( Sd ) ) / 60. )
822	ENDIF
823	IF (minSinceSimStart == 0. .AND. grid%id == 1) THEN
824	PRINT *,' WRF+LMDZ: simulation is starting or comes from a restart!'
825	PRINT *,' since ', minSinceSimStart,' minutes has passed since it started'
826	!! wrftestrst = .TRUE.
827	wrftestrst = .TRUE.
828	ELSE
829	wrftestrst = .FALSE.
830	END IF
831	PRINT *,' Lluis: minSinceSimStart: ', minSinceSimStart
832	! Checking for time for input from auxiliar input 4
833	IF( WRFU_AlarmIsRinging( grid%alarms( first_auxinput + 3 ), rc=rc ) ) THEN
834	PRINT *,' WRF lowbdy time-step!!!!!'
835	wrftestin = .TRUE.
836	CALL WRFU_AlarmRingerOff( grid%alarms( first_auxinput + 3 ), rc=rc )
837	ELSE
838	wrftestin = .FALSE.
839	END IF
840
841	PRINT *,' grid id: ',grid%id
842
843	CALL nl_get_mminlu ( grid%id, mminlu )
844	IF (config_flags%sf_surface_physics == RUCLSMSCHEME) THEN
845	mminsl = 'STAS-RUC'
846	ELSE
847	mminsl = 'STAS'
848	END IF
849
850	grid%qv_2 = moist(:,:,:,P_QV)
851
852	! L. Fita, LMD. July 2014. Getting hour of the day
853	CALL domain_clock_get( grid, current_time=CurrTime )
854	CALL WRFU_TimeGet( CurrTime, YY= julyr, dayOfYear=julday, H=hr, M=minute, S=sec, MS=ms, rc=rc)
855	! Julian day hour (0, 1) !!
856	gmt=(hr+real(minute)/60.+real(sec)/3600.+real(ms)/(1000*3600))/24.
857
858	! Checking for NaNs (should not be necessary but....)
859	im2 = ims + (ime - ims) / 2
860	jm2 = jms + (jme - jms) / 2
861	DO iz = kms, kme
862
863	IF (grid%t_2(im2,iz,jm2) /= grid%t_2(im2,iz,jm2) .OR. ABS(grid%t_2(im2,iz,jm2)) > 10000. ) THEN
864	PRINT *,TRIM(errmsg)
865	WRITE(wrF_err_message,*)'solve_em: wrong T value=', &
866	grid%t_2(im2,iz,jm2),' at: ', im2,', ', iz,', ', jm2,' !!!'
867	#ifdef DM_PARALLEL
868	CALL wrf_error_fatal(TRIM(wrf_err_message))
869	#else
870	PRINT *,TRIM(wrf_err_message)
871	STOP
872	#endif
873	END IF
874	END DO
875	! Checking for NaNs (should not be necessary but....)
876	IF (grid%psfc(im2,jm2) /= grid%psfc(im2,jm2) .OR. ABS(grid%psfc(im2,jm2)) > 1000000. ) THEN
877	PRINT *,errmsg
878	WRITE(wrF_err_message,*)'solve_em: wrong PSFC value=', &
879	grid%psfc(im2,jm2),' at: ', im2 ,', ', jm2, ' !!!'
880	#ifdef DM_PARALLEL
881	CALL wrf_error_fatal(TRIM(wrf_err_message))
882	#else
883	PRINT *,TRIM(wrf_err_message)
884	STOP
885	#endif
886	END IF
887
888	IF (config_flags%lmdz_physics) THEN
889
890	CALL call_lmdz_phys( &
891	& WRF_GRID=grid, WRF_XTIME=grid%xtime, &
892	& WRF_RESTART_ALARM=WRFU_AlarmIsRinging( grid%alarms( RESTART_ALARM ), rc=rc)&
893	& ,WRF_LON = grid%xlong, WRF_LAT=grid%xlat, &
894	& WRF_T=grid%t_2, WRF_U=grid%u_2, WRF_V=grid%v_2, &
895	& WRF_PERP=grid%P, WRF_BASEP=grid%PB, &
896	& WBDYW=config_flags%spec_bdy_width, WRF_ISRESTART=wrftestrst, &
897	& WRF_ISLOWBDYIN=wrftestin &
898	! Dimension arguments
899	& ,WIDS=ids,WIDE=ide, WJDS=jds,WJDE=jde, WKDS=kds,WKDE=kde &
900	& ,WIMS=ims,WIME=ime, WJMS=jms,WJME=jme, WKMS=kms,WKME=kme &
901	& ,WIPS=ips,WIPE=ipe, WJPS=jps,WJPE=jpe, WKPS=kps,WKPE=kpe &
902	& ,WI_START=grid%i_start,WI_END=MIN(grid%i_end, ide-1) &
903	& ,WJ_START=grid%j_start,WJ_END=MIN(grid%j_end, jde-1) &
904	& ,WKTS=k_start, WKTE=MIN(k_end,kde-1) &
905	& ,WNUM_TILES=grid%num_tiles &
906	& ,WNUM3DM=num_3d_m, WPARFIRSTSCAL=PARAM_FIRST_SCALAR, &
907	& WNX=config_flags%e_we, WNY=config_flags%e_sn, &
908	& WNZ=config_flags%e_vert, WJULDAY=FLOAT(julday), WGMT=gmt, &
909	& WTIME_STEP=REAL(config_flags%time_step), &
910	& WRF_FULLETA=grid%znw, WRF_HALFETA=grid%znu, WRF_DFULLETA=grid%dnw, &
911	& WRF_FULLPRES=p8w, WRF_PERGEOPOT=grid%ph_2, &
912	& WRF_BASEGEOPOT=grid%phb, &
913	& WRF_MOIST=grid%moist, WRF_W=grid%w_2, &
914	& WRF_PTOP=config_flags%p_top_requested, &
915	& WRF_PERMASS=grid%mu_1, WRF_BASEMASS=grid%mub, &
916	& WRF_MUT=grid%mut, WRF_MUU=grid%muu, WRF_MUV=grid%muv, &
917	!! & WRF_UTEND=grid%ru_tend, WRF_VTEND=grid%rv_tend, &
918	!! & WRF_TTEND=t_tend, &
919	& WRF_UTEND=ru_tendf, WRF_VTEND=rv_tendf, &
920	& WRF_TTEND=t_tendf, &
921	& WRF_MOISTTEND=moist_tend, WRF_PSFCTEND=grid%dpsdt, &
922	& WRF_QVID=P_QV, WRF_QCID=P_QC, WRF_QRID=P_QR, &
923	& WRF_QSID=P_QS, WRF_QIID=P_QI, WRF_QHID=P_QH, WRF_QGID=P_QG, &
924	! L. Fita. July 2013. Now defined as local dummy variables
925	! & WRF_DUDYN=???????????????, WRF_PVTHETA=????????????????????????, &
926	! & WRF_CLESPHY=?????????????, WRF_PRESNIVS=???????????????????????, &
927	& WRF_MAPFT=grid%msft, WRF_MAPFU=grid%msfu, WRF_MAPFV=grid%msfv, &
928	& WRF_DX=grid%dx, WRF_DY=grid%dy, &
929	& WRF_DBG=model_config_rec%debug_level, LANDCAT=mminlu, &
930	& SOILCAT=mminsl, &
931	& WRF_L_PBL=config_flags%lmdz_iflag_pbl, &
932	& WRF_L_CON=config_flags%lmdz_iflag_con, &
933	& WRF_L_THERMALS=config_flags%lmdz_iflag_thermals, &
934	& WRF_L_WAKE=config_flags%lmdz_iflag_wake, &
935	& wrf_nsoillayers=config_flags%num_soil_layers, &
936	& ICHECK_P=config_flags%i_check_point, &
937	& JCHECK_P=config_flags%j_check_point, &
938	& KCHECK_P=config_flags%k_check_point &
939	& )
940	END IF
941	#endif
942
943	#ifdef LMDZ1
944	WRITE(message, *)' dyn_em: post lmdz pre step2 rk_step:', rk_step
945	CALL wrf_debug(200, message)
946	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
947	' u_tend: ', ru_tendf(im2,1,jm2)
948	CALL wrf_debug(200, message)
949	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
950	'p sfc: ',p8w(im2,kms,jm2)
951	CALL wrf_debug(200, message)
952	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
953	CALL wrf_debug(200, message)
954	#endif
955
956	CALL first_rk_step_part2 ( grid, config_flags &
957	, moist , moist_tend &
958	, chem , chem_tend &
959	, tracer, tracer_tend &
960	, scalar , scalar_tend &
961	, fdda3d, fdda2d &
962	, ru_tendf, rv_tendf &
963	, rw_tendf, t_tendf &
964	, ph_tendf, mu_tendf &
965	, tke_tend &
966	, adapt_step_flag , curr_secs &
967	, psim , psih , wspd , gz1oz0 &
968	, br , chklowq &
969	, cu_act_flag , hol , th_phy &
970	, pi_phy , p_phy , grid%t_phy &
971	, u_phy , v_phy &
972	, dz8w , p8w , t8w , rho_phy , rho &
973	, nba_mij, num_nba_mij & !JDM
974	, nba_rij, num_nba_rij & !JDM
975	, ids, ide, jds, jde, kds, kde &
976	, ims, ime, jms, jme, kms, kme &
977	, ips, ipe, jps, jpe, kps, kpe &
978	, imsx, imex, jmsx, jmex, kmsx, kmex &
979	, ipsx, ipex, jpsx, jpex, kpsx, kpex &
980	, imsy, imey, jmsy, jmey, kmsy, kmey &
981	, ipsy, ipey, jpsy, jpey, kpsy, kpey &
982	, k_start , k_end &
983	)
984
985	END IF rk_step_is_one
986
987	#ifdef LMDZ1
988	WRITE(message, *)' dyn_em: post step2 pre rk_tendency rk_step: ',rk_step
989	CALL wrf_debug(200, message)
990	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
991	' u_tend: ', ru_tendf(im2,1,jm2)
992	CALL wrf_debug(200, message)
993	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
994	'p sfc: ',p8w(im2,kms,jm2)
995	CALL wrf_debug(200, message)
996	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
997	CALL wrf_debug(200, message)
998	#endif
999
1000	BENCH_START(rk_tend_tim)
1001	!$OMP PARALLEL DO &
1002	!$OMP PRIVATE ( ij )
1003	DO ij = 1 , grid%num_tiles
1004
1005	CALL wrf_debug ( 200 , ' call rk_tendency' )
1006	CALL rk_tendency ( config_flags, rk_step &
1007	,grid%ru_tend, grid%rv_tend, rw_tend, ph_tend, t_tend &
1008	,ru_tendf, rv_tendf, rw_tendf, ph_tendf, t_tendf &
1009	,mu_tend, grid%u_save, grid%v_save, w_save, ph_save &
1010	,grid%t_save, mu_save, grid%rthften &
1011	,grid%ru, grid%rv, grid%rw, grid%ww &
1012	,grid%u_2, grid%v_2, grid%w_2, grid%t_2, grid%ph_2 &
1013	,grid%u_1, grid%v_1, grid%w_1, grid%t_1, grid%ph_1 &
1014	,grid%h_diabatic, grid%phb, grid%t_init &
1015	,grid%mu_2, grid%mut, grid%muu, grid%muv, grid%mub &
1016	,grid%al, grid%alt, grid%p, grid%pb, grid%php, cqu, cqv, cqw &
1017	,grid%u_base, grid%v_base, grid%t_base, grid%qv_base, grid%z_base &
1018	,grid%msfux,grid%msfuy, grid%msfvx, grid%msfvx_inv &
1019	,grid%msfvy, grid%msftx,grid%msfty, grid%clat, grid%f, grid%e, grid%sina, grid%cosa &
1020	,grid%fnm, grid%fnp, grid%rdn, grid%rdnw &
1021	,grid%dt, grid%rdx, grid%rdy, grid%khdif, grid%kvdif, grid%xkmh, grid%xkhh &
1022	,grid%diff_6th_opt, grid%diff_6th_factor &
1023	,grid%dampcoef,grid%zdamp,config_flags%damp_opt,config_flags%rad_nudge &
1024	,grid%cf1, grid%cf2, grid%cf3, grid%cfn, grid%cfn1, num_3d_m &
1025	,config_flags%non_hydrostatic, config_flags%top_lid &
1026	,grid%u_frame, grid%v_frame &
1027	,ids, ide, jds, jde, kds, kde &
1028	,ims, ime, jms, jme, kms, kme &
1029	,grid%i_start(ij), grid%i_end(ij) &
1030	,grid%j_start(ij), grid%j_end(ij) &
1031	,k_start, k_end &
1032	,max_vert_cfl_tmp(ij), max_horiz_cfl_tmp(ij) )
1033
1034	END DO
1035	#ifdef LMDZ1
1036	WRITE(message, *)' dyn_em: post rk_tendency'
1037	CALL wrf_debug(200, message)
1038	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1039	' u_tend: ', ru_tendf(im2,1,jm2)
1040	CALL wrf_debug(200, message)
1041	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1042	'p sfc: ',p8w(im2,kms,jm2)
1043	CALL wrf_debug(200, message)
1044	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1045	CALL wrf_debug(200, message)
1046	#endif
1047	!$OMP END PARALLEL DO
1048	BENCH_END(rk_tend_tim)
1049
1050	IF (config_flags%use_adaptive_time_step) THEN
1051	DO ij = 1 , grid%num_tiles
1052	IF (max_horiz_cfl_tmp(ij) .GT. grid%max_horiz_cfl) THEN
1053	grid%max_horiz_cfl = max_horiz_cfl_tmp(ij)
1054	ENDIF
1055	IF (max_vert_cfl_tmp(ij) .GT. grid%max_vert_cfl) THEN
1056	grid%max_vert_cfl = max_vert_cfl_tmp(ij)
1057	ENDIF
1058	END DO
1059
1060	IF (grid%max_horiz_cfl .GT. grid%max_cfl_val) THEN
1061	grid%max_cfl_val = grid%max_horiz_cfl
1062	ENDIF
1063	IF (grid%max_vert_cfl .GT. grid%max_cfl_val) THEN
1064	grid%max_cfl_val = grid%max_vert_cfl
1065	ENDIF
1066	ENDIF
1067
1068	BENCH_START(relax_bdy_dry_tim)
1069	!$OMP PARALLEL DO &
1070	!$OMP PRIVATE ( ij )
1071	DO ij = 1 , grid%num_tiles
1072
1073	IF( (config_flags%specified .or. config_flags%nested) .and. rk_step == 1 ) THEN
1074
1075	CALL relax_bdy_dry ( config_flags, &
1076	grid%u_save, grid%v_save, ph_save, grid%t_save, &
1077	w_save, mu_tend, &
1078	grid%ru, grid%rv, grid%ph_2, grid%t_2, &
1079	grid%w_2, grid%mu_2, grid%mut, &
1080	grid%u_bxs,grid%u_bxe,grid%u_bys,grid%u_bye, &
1081	grid%v_bxs,grid%v_bxe,grid%v_bys,grid%v_bye, &
1082	grid%ph_bxs,grid%ph_bxe,grid%ph_bys,grid%ph_bye, &
1083	grid%t_bxs,grid%t_bxe,grid%t_bys,grid%t_bye, &
1084	grid%w_bxs,grid%w_bxe,grid%w_bys,grid%w_bye, &
1085	grid%mu_bxs,grid%mu_bxe,grid%mu_bys,grid%mu_bye, &
1086	grid%u_btxs,grid%u_btxe,grid%u_btys,grid%u_btye, &
1087	grid%v_btxs,grid%v_btxe,grid%v_btys,grid%v_btye, &
1088	grid%ph_btxs,grid%ph_btxe,grid%ph_btys,grid%ph_btye, &
1089	grid%t_btxs,grid%t_btxe,grid%t_btys,grid%t_btye, &
1090	grid%w_btxs,grid%w_btxe,grid%w_btys,grid%w_btye, &
1091	grid%mu_btxs,grid%mu_btxe,grid%mu_btys,grid%mu_btye, &
1092	config_flags%spec_bdy_width, grid%spec_zone, grid%relax_zone, &
1093	grid%dtbc, grid%fcx, grid%gcx, &
1094	ids,ide, jds,jde, kds,kde, &
1095	ims,ime, jms,jme, kms,kme, &
1096	ips,ipe, jps,jpe, kps,kpe, &
1097	grid%i_start(ij), grid%i_end(ij), &
1098	grid%j_start(ij), grid%j_end(ij), &
1099	k_start, k_end )
1100
1101	ENDIF
1102
1103	CALL rk_addtend_dry( grid%ru_tend, grid%rv_tend, rw_tend, ph_tend, t_tend, &
1104	ru_tendf, rv_tendf, rw_tendf, ph_tendf, t_tendf, &
1105	grid%u_save, grid%v_save, w_save, ph_save, grid%t_save, &
1106	mu_tend, mu_tendf, rk_step, &
1107	grid%h_diabatic, grid%mut, grid%msftx, &
1108	grid%msfty, grid%msfux,grid%msfuy, &
1109	grid%msfvx, grid%msfvx_inv, grid%msfvy, &
1110	ids,ide, jds,jde, kds,kde, &
1111	ims,ime, jms,jme, kms,kme, &
1112	ips,ipe, jps,jpe, kps,kpe, &
1113	grid%i_start(ij), grid%i_end(ij), &
1114	grid%j_start(ij), grid%j_end(ij), &
1115	k_start, k_end )
1116
1117	IF( config_flags%specified .or. config_flags%nested ) THEN
1118	CALL spec_bdy_dry ( config_flags, &
1119	grid%ru_tend, grid%rv_tend, ph_tend, t_tend, &
1120	rw_tend, mu_tend, &
1121	grid%u_bxs,grid%u_bxe,grid%u_bys,grid%u_bye, &
1122	grid%v_bxs,grid%v_bxe,grid%v_bys,grid%v_bye, &
1123	grid%ph_bxs,grid%ph_bxe,grid%ph_bys,grid%ph_bye, &
1124	grid%t_bxs,grid%t_bxe,grid%t_bys,grid%t_bye, &
1125	grid%w_bxs,grid%w_bxe,grid%w_bys,grid%w_bye, &
1126	grid%mu_bxs,grid%mu_bxe,grid%mu_bys,grid%mu_bye, &
1127	grid%u_btxs,grid%u_btxe,grid%u_btys,grid%u_btye, &
1128	grid%v_btxs,grid%v_btxe,grid%v_btys,grid%v_btye, &
1129	grid%ph_btxs,grid%ph_btxe,grid%ph_btys,grid%ph_btye, &
1130	grid%t_btxs,grid%t_btxe,grid%t_btys,grid%t_btye, &
1131	grid%w_btxs,grid%w_btxe,grid%w_btys,grid%w_btye, &
1132	grid%mu_btxs,grid%mu_btxe,grid%mu_btys,grid%mu_btye, &
1133	config_flags%spec_bdy_width, grid%spec_zone, &
1134	ids,ide, jds,jde, kds,kde, & ! domain dims
1135	ims,ime, jms,jme, kms,kme, & ! memory dims
1136	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
1137	grid%i_start(ij), grid%i_end(ij), &
1138	grid%j_start(ij), grid%j_end(ij), &
1139	k_start, k_end )
1140
1141	ENDIF
1142
1143	END DO
1144	!$OMP END PARALLEL DO
1145	BENCH_END(relax_bdy_dry_tim)
1146
1147	!<DESCRIPTION>
1148	!<pre>
1149	! (3) Small (acoustic,sound) steps.
1150	!
1151	! Several acoustic steps are taken each RK pass. A small step
1152	! sequence begins with calculating perturbation variables
1153	! and coupling them to the column dry-air-mass mu
1154	! (call to small_step_prep). This is followed by computing
1155	! coefficients for the vertically implicit part of the
1156	! small timestep (call to calc_coef_w).
1157	!
1158	! The small steps are taken
1159	! in the named loop "small_steps:". In the small_steps loop, first
1160	! the horizontal momentum (u and v) are advanced (call to advance_uv),
1161	! next mu and theta are advanced (call to advance_mu_t) followed by
1162	! advancing w and the geopotential (call to advance_w). Diagnostic
1163	! values for pressure and inverse density are updated at the end of
1164	! each small_step.
1165	!
1166	! The small-step section ends with the change of the perturbation variables
1167	! back to full variables (call to small_step_finish).
1168	!</pre>
1169	!</DESCRIPTION>
1170
1171	BENCH_START(small_step_prep_tim)
1172	!$OMP PARALLEL DO &
1173	!$OMP PRIVATE ( ij )
1174	DO ij = 1 , grid%num_tiles
1175
1176	! Calculate coefficients for the vertically implicit acoustic/gravity wave
1177	! integration. We only need calculate these for the first pass through -
1178	! the predictor step. They are reused as is for the corrector step.
1179	! For third-order RK, we need to recompute these after the first
1180	! predictor because we may have changed the small timestep -> grid%dts.
1181
1182	CALL wrf_debug ( 200 , ' call small_step_prep ' )
1183
1184	#ifdef LMDZ1
1185	WRITE(message, *)' dyn_em: before small_step_prep'
1186	CALL wrf_debug(200, message)
1187	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1188	' u_tend: ', ru_tendf(im2,1,jm2)
1189	CALL wrf_debug(200, message)
1190	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1191	'p sfc: ',p8w(im2,kms,jm2)
1192	CALL wrf_debug(200, message)
1193	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1194	CALL wrf_debug(200, message)
1195	#endif
1196	CALL small_step_prep( grid%u_1,grid%u_2,grid%v_1,grid%v_2,grid%w_1,grid%w_2, &
1197	grid%t_1,grid%t_2,grid%ph_1,grid%ph_2, &
1198	grid%mub, grid%mu_1, grid%mu_2, &
1199	grid%muu, muus, grid%muv, muvs, &
1200	grid%mut, grid%muts, grid%mudf, &
1201	grid%u_save, grid%v_save, w_save, &
1202	grid%t_save, ph_save, mu_save, &
1203	grid%ww, ww1, &
1204	grid%dnw, c2a, grid%pb, grid%p, grid%alt, &
1205	grid%msfux,grid%msfuy, grid%msfvx, grid%msfvx_inv, &
1206	grid%msfvy, grid%msftx,grid%msfty, &
1207	grid%rdx, grid%rdy, rk_step, &
1208	ids, ide, jds, jde, kds, kde, &
1209	ims, ime, jms, jme, kms, kme, &
1210	grid%i_start(ij), grid%i_end(ij), &
1211	grid%j_start(ij), grid%j_end(ij), &
1212	k_start , k_end )
1213
1214	#ifdef LMDZ1
1215	WRITE(message, *)' dyn_em: post small_step_prep'
1216	CALL wrf_debug(200, message)
1217	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1218	' u_tend: ', ru_tendf(im2,1,jm2)
1219	CALL wrf_debug(200, message)
1220	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1221	'p sfc: ',p8w(im2,kms,jm2)
1222	CALL wrf_debug(200, message)
1223	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1224	CALL wrf_debug(200, message)
1225	#endif
1226	CALL calc_p_rho( grid%al, grid%p, grid%ph_2, &
1227	grid%alt, grid%t_2, grid%t_save, c2a, pm1, &
1228	grid%mu_2, grid%muts, grid%znu, t0, &
1229	grid%rdnw, grid%dnw, grid%smdiv, &
1230	config_flags%non_hydrostatic, 0, &
1231	ids, ide, jds, jde, kds, kde, &
1232	ims, ime, jms, jme, kms, kme, &
1233	grid%i_start(ij), grid%i_end(ij), &
1234	grid%j_start(ij), grid%j_end(ij), &
1235	k_start , k_end )
1236	#ifdef LMDZ1
1237	WRITE(message, *)' dyn_em: post calc_p_rho'
1238	CALL wrf_debug(200, message)
1239	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1240	' u_tend: ', ru_tendf(im2,1,jm2)
1241	CALL wrf_debug(200, message)
1242	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1243	'p sfc: ',p8w(im2,kms,jm2)
1244	CALL wrf_debug(200, message)
1245	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1246	CALL wrf_debug(200, message)
1247	#endif
1248
1249	IF (config_flags%non_hydrostatic) THEN
1250	CALL calc_coef_w( a,alpha,gamma, &
1251	grid%mut, cqw, &
1252	grid%rdn, grid%rdnw, c2a, &
1253	dts_rk, g, grid%epssm, &
1254	config_flags%top_lid, &
1255	ids, ide, jds, jde, kds, kde, &
1256	ims, ime, jms, jme, kms, kme, &
1257	grid%i_start(ij), grid%i_end(ij), &
1258	grid%j_start(ij), grid%j_end(ij), &
1259	k_start , k_end )
1260	ENDIF
1261	#ifdef LMDZ1
1262	WRITE(message, *)' dyn_em: post calc_coef_w'
1263	CALL wrf_debug(200, message)
1264	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1265	' u_tend: ', ru_tendf(im2,1,jm2)
1266	CALL wrf_debug(200, message)
1267	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1268	'p sfc: ',p8w(im2,kms,jm2)
1269	CALL wrf_debug(200, message)
1270	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1271	CALL wrf_debug(200, message)
1272	WRITE(message,*)' al: ',grid%al(im2,km2,jm2), ' p: ', grid%p(im2,1,jm2), &
1273	' ph: ',grid%ph_2(im2,1,jm2)
1274	CALL wrf_debug(200, message)
1275	WRITE(message,*)' mu: ',grid%mu_2(im2,jm2), 'alt: ', grid%alt(im2,1,jm2), &
1276	' mu: ',grid%mu_2(im2,jm2), ' znu: ', grid%znu(1),' ph 1: ',grid%ph_2(im2,2,jm2)
1277	CALL wrf_debug(200, message)
1278	PRINT *,' c2a: ',c2a(im2,1,jm2), 't_2: ', grid%t_2(im2,1,jm2), &
1279	' t_save: ',grid%t_save(im2,1,jm2),' pm1: ',pm1(im2,1,jm2)
1280	CALL wrf_debug(200, message)
1281	#endif
1282	ENDDO
1283	!$OMP END PARALLEL DO
1284	BENCH_END(small_step_prep_tim)
1285
1286	#ifdef DM_PARALLEL
1287	!-----------------------------------------------------------------------
1288	! Stencils for patch communications (WCS, 29 June 2001)
1289	! Note: the small size of this halo exchange reflects the
1290	! fact that we are carrying the uncoupled variables
1291	! as state variables in the mass coordinate model, as
1292	! opposed to the coupled variables as in the height
1293	! coordinate model.
1294	!
1295	! * * * * *
1296	! * * * * * * * * *
1297	! * + * * + * * * + * *
1298	! * * * * * * * * *
1299	! * * * * *
1300	!
1301	! 3D variables - note staggering! ph_2(Z), u_save(X), v_save(Y)
1302	!
1303	! ph_2 x
1304	! al x
1305	! p x
1306	! t_1 x
1307	! t_save x
1308	! u_save x
1309	! v_save x
1310	!
1311	! the following are 2D (xy) variables
1312	!
1313	! mu_1 x
1314	! mu_2 x
1315	! mudf x
1316	! php x
1317	! alt x
1318	! pb x
1319	!--------------------------------------------------------------
1320	# include "HALO_EM_B.inc"
1321	# include "PERIOD_BDY_EM_B.inc"
1322	#endif
1323
1324	BENCH_START(set_phys_bc2_tim)
1325	!$OMP PARALLEL DO &
1326	!$OMP PRIVATE ( ij )
1327
1328	DO ij = 1 , grid%num_tiles
1329
1330	CALL set_physical_bc3d( grid%ru_tend, 'u', config_flags, &
1331	ids, ide, jds, jde, kds, kde, &
1332	ims, ime, jms, jme, kms, kme, &
1333	ips, ipe, jps, jpe, kps, kpe, &
1334	grid%i_start(ij), grid%i_end(ij), &
1335	grid%j_start(ij), grid%j_end(ij), &
1336	k_start , k_end )
1337
1338	CALL set_physical_bc3d( grid%rv_tend, 'v', config_flags, &
1339	ids, ide, jds, jde, kds, kde, &
1340	ims, ime, jms, jme, kms, kme, &
1341	ips, ipe, jps, jpe, kps, kpe, &
1342	grid%i_start(ij), grid%i_end(ij), &
1343	grid%j_start(ij), grid%j_end(ij), &
1344	k_start , k_end )
1345
1346	CALL set_physical_bc3d( grid%ph_2, 'w', config_flags, &
1347	ids, ide, jds, jde, kds, kde, &
1348	ims, ime, jms, jme, kms, kme, &
1349	ips, ipe, jps, jpe, kps, kpe, &
1350	grid%i_start(ij), grid%i_end(ij), &
1351	grid%j_start(ij), grid%j_end(ij), &
1352	k_start , k_end )
1353
1354	CALL set_physical_bc3d( grid%al, 'p', config_flags, &
1355	ids, ide, jds, jde, kds, kde, &
1356	ims, ime, jms, jme, kms, kme, &
1357	ips, ipe, jps, jpe, kps, kpe, &
1358	grid%i_start(ij), grid%i_end(ij), &
1359	grid%j_start(ij), grid%j_end(ij), &
1360	k_start , k_end )
1361
1362	CALL set_physical_bc3d( grid%p, 'p', config_flags, &
1363	ids, ide, jds, jde, kds, kde, &
1364	ims, ime, jms, jme, kms, kme, &
1365	ips, ipe, jps, jpe, kps, kpe, &
1366	grid%i_start(ij), grid%i_end(ij), &
1367	grid%j_start(ij), grid%j_end(ij), &
1368	k_start , k_end )
1369
1370	CALL set_physical_bc3d( grid%t_1, 'p', config_flags, &
1371	ids, ide, jds, jde, kds, kde, &
1372	ims, ime, jms, jme, kms, kme, &
1373	ips, ipe, jps, jpe, kps, kpe, &
1374	grid%i_start(ij), grid%i_end(ij), &
1375	grid%j_start(ij), grid%j_end(ij), &
1376	k_start , k_end )
1377
1378	CALL set_physical_bc3d( grid%t_save, 't', config_flags, &
1379	ids, ide, jds, jde, kds, kde, &
1380	ims, ime, jms, jme, kms, kme, &
1381	ips, ipe, jps, jpe, kps, kpe, &
1382	grid%i_start(ij), grid%i_end(ij), &
1383	grid%j_start(ij), grid%j_end(ij), &
1384	k_start , k_end )
1385
1386	CALL set_physical_bc2d( grid%mu_1, 't', config_flags, &
1387	ids, ide, jds, jde, &
1388	ims, ime, jms, jme, &
1389	ips, ipe, jps, jpe, &
1390	grid%i_start(ij), grid%i_end(ij), &
1391	grid%j_start(ij), grid%j_end(ij) )
1392
1393	CALL set_physical_bc2d( grid%mu_2, 't', config_flags, &
1394	ids, ide, jds, jde, &
1395	ims, ime, jms, jme, &
1396	ips, ipe, jps, jpe, &
1397	grid%i_start(ij), grid%i_end(ij), &
1398	grid%j_start(ij), grid%j_end(ij) )
1399
1400	CALL set_physical_bc2d( grid%mudf, 't', config_flags, &
1401	ids, ide, jds, jde, &
1402	ims, ime, jms, jme, &
1403	ips, ipe, jps, jpe, &
1404	grid%i_start(ij), grid%i_end(ij), &
1405	grid%j_start(ij), grid%j_end(ij) )
1406
1407	END DO
1408	!$OMP END PARALLEL DO
1409	BENCH_END(set_phys_bc2_tim)
1410	small_steps : DO iteration = 1 , number_of_small_timesteps
1411
1412	! Boundary condition time (or communication time).
1413	#ifdef DM_PARALLEL
1414	# include "PERIOD_BDY_EM_B.inc"
1415	#endif
1416
1417	!$OMP PARALLEL DO &
1418	!$OMP PRIVATE ( ij )
1419
1420	DO ij = 1 , grid%num_tiles
1421
1422	BENCH_START(advance_uv_tim)
1423	#ifdef LMDZ1
1424	WRITE(message, *)' dyn_em: before advance_uv'
1425	CALL wrf_debug(200, message)
1426	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1427	' u_tend: ', ru_tendf(im2,1,jm2)
1428	CALL wrf_debug(200, message)
1429	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1430	'p sfc: ',p8w(im2,kms,jm2)
1431	CALL wrf_debug(200, message)
1432	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1433	CALL wrf_debug(200, message)
1434	WRITE(message,*)' al: ',grid%al(im2,km2,jm2), ' p: ', grid%p(im2,1,jm2), &
1435	' ph: ',grid%ph_2(im2,1,jm2)
1436	CALL wrf_debug(200, message)
1437	WRITE(message,*)' mu: ',grid%mu_2(im2,jm2), 'alt: ', grid%alt(im2,1,jm2), &
1438	' mu: ',grid%mu_2(im2,jm2), ' znu: ', grid%znu(1),' ph 1: ',grid%ph_2(im2,2,jm2)
1439	CALL wrf_debug(200, message)
1440	PRINT *,' c2a: ',c2a(im2,1,jm2), 't_2: ', grid%t_2(im2,1,jm2), &
1441	' t_save: ',grid%t_save(im2,1,jm2),' pm1: ',pm1(im2,1,jm2)
1442	CALL wrf_debug(200, message)
1443	#endif
1444
1445	CALL advance_uv ( grid%u_2, grid%ru_tend, grid%v_2, grid%rv_tend, &
1446	grid%p, grid%pb, &
1447	grid%ph_2, grid%php, grid%alt, grid%al, &
1448	grid%mu_2, &
1449	grid%muu, cqu, grid%muv, cqv, grid%mudf, &
1450	grid%msfux, grid%msfuy, grid%msfvx, &
1451	grid%msfvx_inv, grid%msfvy, &
1452	grid%rdx, grid%rdy, dts_rk, &
1453	grid%cf1, grid%cf2, grid%cf3, grid%fnm, grid%fnp, &
1454	grid%emdiv, &
1455	grid%rdnw, config_flags,grid%spec_zone, &
1456	config_flags%non_hydrostatic, config_flags%top_lid, &
1457	ids, ide, jds, jde, kds, kde, &
1458	ims, ime, jms, jme, kms, kme, &
1459	grid%i_start(ij), grid%i_end(ij), &
1460	grid%j_start(ij), grid%j_end(ij), &
1461	k_start , k_end )
1462	BENCH_END(advance_uv_tim)
1463
1464	END DO
1465	!$OMP END PARALLEL DO
1466	#ifdef LMDZ1
1467	WRITE(message, *)' dyn_em: after advance_uv'
1468	CALL wrf_debug(200, message)
1469	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1470	' u_tend: ', ru_tendf(im2,1,jm2)
1471	CALL wrf_debug(200, message)
1472	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1473	'p sfc: ',p8w(im2,kms,jm2)
1474	CALL wrf_debug(200, message)
1475	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1476	CALL wrf_debug(200, message)
1477	WRITE(message,*)' al: ',grid%al(im2,km2,jm2), ' p: ', grid%p(im2,1,jm2), &
1478	' ph: ',grid%ph_2(im2,1,jm2)
1479	CALL wrf_debug(200, message)
1480	WRITE(message,*)' mu: ',grid%mu_2(im2,jm2), 'alt: ', grid%alt(im2,1,jm2), &
1481	' mu: ',grid%mu_2(im2,jm2), ' znu: ', grid%znu(1),' ph 1: ',grid%ph_2(im2,2,jm2)
1482	CALL wrf_debug(200, message)
1483	PRINT *,' c2a: ',c2a(im2,1,jm2), 't_2: ', grid%t_2(im2,1,jm2), &
1484	' t_save: ',grid%t_save(im2,1,jm2),' pm1: ',pm1(im2,1,jm2)
1485	CALL wrf_debug(200, message)
1486	#endif
1487
1488	!-----------------------------------------------------------
1489	! acoustic integration polar filter for smallstep u, v
1490	!-----------------------------------------------------------
1491
1492	IF (config_flags%polar) THEN
1493
1494	CALL pxft ( grid=grid &
1495	,lineno=__LINE__ &
1496	,flag_uv = 1 &
1497	,flag_rurv = 0 &
1498	,flag_wph = 0 &
1499	,flag_ww = 0 &
1500	,flag_t = 0 &
1501	,flag_mu = 0 &
1502	,flag_mut = 0 &
1503	,flag_moist = 0 &
1504	,flag_chem = 0 &
1505	,flag_tracer = 0 &
1506	,flag_scalar = 0 &
1507	,positive_definite = .FALSE. &
1508	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
1509	,fft_filter_lat = config_flags%fft_filter_lat &
1510	,dclat = dclat &
1511	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
1512	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
1513	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
1514	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
1515	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
1516
1517	END IF
1518
1519	!-----------------------------------------------------------
1520	! end acoustic integration polar filter for smallstep u, v
1521	!-----------------------------------------------------------
1522
1523	!$OMP PARALLEL DO &
1524	!$OMP PRIVATE ( ij )
1525	DO ij = 1 , grid%num_tiles
1526
1527	BENCH_START(spec_bdy_uv_tim)
1528	IF( config_flags%specified .or. config_flags%nested ) THEN
1529	CALL spec_bdyupdate(grid%u_2, grid%ru_tend, dts_rk, &
1530	'u' , config_flags, &
1531	grid%spec_zone, &
1532	ids,ide, jds,jde, kds,kde, & ! domain dims
1533	ims,ime, jms,jme, kms,kme, & ! memory dims
1534	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
1535	grid%i_start(ij), grid%i_end(ij), &
1536	grid%j_start(ij), grid%j_end(ij), &
1537	k_start , k_end )
1538
1539	CALL spec_bdyupdate(grid%v_2, grid%rv_tend, dts_rk, &
1540	'v' , config_flags, &
1541	grid%spec_zone, &
1542	ids,ide, jds,jde, kds,kde, & ! domain dims
1543	ims,ime, jms,jme, kms,kme, & ! memory dims
1544	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
1545	grid%i_start(ij), grid%i_end(ij), &
1546	grid%j_start(ij), grid%j_end(ij), &
1547	k_start , k_end )
1548
1549	ENDIF
1550	BENCH_END(spec_bdy_uv_tim)
1551
1552	END DO
1553	!$OMP END PARALLEL DO
1554
1555	#ifdef DM_PARALLEL
1556	!
1557	! Stencils for patch communications (WCS, 29 June 2001)
1558	!
1559	! * *
1560	! * + * * + * +
1561	! * *
1562	!
1563	! u_2 x
1564	! v_2 x
1565	!
1566	# include "HALO_EM_C.inc"
1567	#endif
1568
1569	!$OMP PARALLEL DO &
1570	!$OMP PRIVATE ( ij )
1571	DO ij = 1 , grid%num_tiles
1572
1573	! advance the mass in the column, theta, and calculate ww
1574
1575	BENCH_START(advance_mu_t_tim)
1576	CALL advance_mu_t( grid%ww, ww1, grid%u_2, grid%u_save, grid%v_2, grid%v_save, &
1577	grid%mu_2, grid%mut, muave, grid%muts, grid%muu, grid%muv, &
1578	grid%mudf, grid%ru_m, grid%rv_m, grid%ww_m, &
1579	grid%t_2, grid%t_save, t_2save, t_tend, &
1580	mu_tend, &
1581	grid%rdx, grid%rdy, dts_rk, grid%epssm, &
1582	grid%dnw, grid%fnm, grid%fnp, grid%rdnw, &
1583	grid%msfux,grid%msfuy, grid%msfvx, grid%msfvx_inv, &
1584	grid%msfvy, grid%msftx,grid%msfty, &
1585	iteration, config_flags, &
1586	ids, ide, jds, jde, kds, kde, &
1587	ims, ime, jms, jme, kms, kme, &
1588	grid%i_start(ij), grid%i_end(ij), &
1589	grid%j_start(ij), grid%j_end(ij), &
1590	k_start , k_end )
1591	BENCH_END(advance_mu_t_tim)
1592	ENDDO
1593	!$OMP END PARALLEL DO
1594	#ifdef LMDZ1
1595	WRITE(message, *)' dyn_em: after advance_mut'
1596	CALL wrf_debug(200, message)
1597	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1598	' u_tend: ', ru_tendf(im2,1,jm2)
1599	CALL wrf_debug(200, message)
1600	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1601	'p sfc: ',p8w(im2,kms,jm2)
1602	CALL wrf_debug(200, message)
1603	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1604	CALL wrf_debug(200, message)
1605	WRITE(message,*)' al: ',grid%al(im2,km2,jm2), ' p: ', grid%p(im2,1,jm2), &
1606	' ph: ',grid%ph_2(im2,1,jm2)
1607	CALL wrf_debug(200, message)
1608	WRITE(message,*)' mu: ',grid%mu_2(im2,jm2), 'alt: ', grid%alt(im2,1,jm2), &
1609	' mu: ',grid%mu_2(im2,jm2), ' znu: ', grid%znu(1),' ph 1: ',grid%ph_2(im2,2,jm2)
1610	CALL wrf_debug(200, message)
1611	PRINT *,' c2a: ',c2a(im2,1,jm2), 't_2: ', grid%t_2(im2,1,jm2), &
1612	' t_save: ',grid%t_save(im2,1,jm2),' pm1: ',pm1(im2,1,jm2)
1613	CALL wrf_debug(200, message)
1614	#endif
1615
1616	!-----------------------------------------------------------
1617	! acoustic integration polar filter for smallstep mu, t
1618	!-----------------------------------------------------------
1619
1620	IF ( (config_flags%polar) ) THEN
1621
1622	CALL pxft ( grid=grid &
1623	,lineno=__LINE__ &
1624	,flag_uv = 0 &
1625	,flag_rurv = 0 &
1626	,flag_wph = 0 &
1627	,flag_ww = 0 &
1628	,flag_t = 1 &
1629	,flag_mu = 1 &
1630	,flag_mut = 0 &
1631	,flag_moist = 0 &
1632	,flag_chem = 0 &
1633	,flag_tracer = 0 &
1634	,flag_scalar = 0 &
1635	,positive_definite = .FALSE. &
1636	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
1637	,fft_filter_lat = config_flags%fft_filter_lat &
1638	,dclat = dclat &
1639	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
1640	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
1641	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
1642	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
1643	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
1644
1645	grid%muts = grid%mut + grid%mu_2 ! reset muts using filtered mu_2
1646
1647	END IF
1648
1649	!-----------------------------------------------------------
1650	! end acoustic integration polar filter for smallstep mu, t
1651	!-----------------------------------------------------------
1652
1653	BENCH_START(spec_bdy_t_tim)
1654
1655	!$OMP PARALLEL DO &
1656	!$OMP PRIVATE ( ij )
1657	DO ij = 1 , grid%num_tiles
1658
1659	IF( config_flags%specified .or. config_flags%nested ) THEN
1660
1661	CALL spec_bdyupdate(grid%t_2, t_tend, dts_rk, &
1662	't' , config_flags, &
1663	grid%spec_zone, &
1664	ids,ide, jds,jde, kds,kde, &
1665	ims,ime, jms,jme, kms,kme, &
1666	ips,ipe, jps,jpe, kps,kpe, &
1667	grid%i_start(ij), grid%i_end(ij),&
1668	grid%j_start(ij), grid%j_end(ij),&
1669	k_start , k_end )
1670
1671	CALL spec_bdyupdate(grid%mu_2, mu_tend, dts_rk, &
1672	'm' , config_flags, &
1673	grid%spec_zone, &
1674	ids,ide, jds,jde, 1 ,1 , &
1675	ims,ime, jms,jme, 1 ,1 , &
1676	ips,ipe, jps,jpe, 1 ,1 , &
1677	grid%i_start(ij), grid%i_end(ij),&
1678	grid%j_start(ij), grid%j_end(ij),&
1679	1 , 1 )
1680
1681	CALL spec_bdyupdate(grid%muts, mu_tend, dts_rk, &
1682	'm' , config_flags, &
1683	grid%spec_zone, &
1684	ids,ide, jds,jde, 1 ,1 , & ! domain dims
1685	ims,ime, jms,jme, 1 ,1 , & ! memory dims
1686	ips,ipe, jps,jpe, 1 ,1 , & ! patch dims
1687	grid%i_start(ij), grid%i_end(ij), &
1688	grid%j_start(ij), grid%j_end(ij), &
1689	1 , 1 )
1690	ENDIF
1691	BENCH_END(spec_bdy_t_tim)
1692
1693	! small (acoustic) step for the vertical momentum,
1694	! density and coupled potential temperature.
1695
1696
1697	BENCH_START(advance_w_tim)
1698	IF ( config_flags%non_hydrostatic ) THEN
1699	CALL advance_w( grid%w_2, rw_tend, grid%ww, w_save, &
1700	grid%u_2, grid%v_2, &
1701	grid%mu_2, grid%mut, muave, grid%muts, &
1702	t_2save, grid%t_2, grid%t_save, &
1703	grid%ph_2, ph_save, grid%phb, ph_tend, &
1704	grid%ht, c2a, cqw, grid%alt, grid%alb, &
1705	a, alpha, gamma, &
1706	grid%rdx, grid%rdy, dts_rk, t0, grid%epssm, &
1707	grid%dnw, grid%fnm, grid%fnp, grid%rdnw, &
1708	grid%rdn, grid%cf1, grid%cf2, grid%cf3, &
1709	grid%msftx, grid%msfty, &
1710	config_flags, config_flags%top_lid, &
1711	ids,ide, jds,jde, kds,kde, &
1712	ims,ime, jms,jme, kms,kme, &
1713	grid%i_start(ij), grid%i_end(ij), &
1714	grid%j_start(ij), grid%j_end(ij), &
1715	k_start , k_end )
1716	ENDIF
1717	BENCH_END(advance_w_tim)
1718
1719	ENDDO
1720	!$OMP END PARALLEL DO
1721	#ifdef LMDZ1
1722	WRITE(message, *)' dyn_em: after advance_w'
1723	CALL wrf_debug(200, message)
1724	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1725	' u_tend: ', ru_tendf(im2,1,jm2)
1726	CALL wrf_debug(200, message)
1727	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1728	'p sfc: ',p8w(im2,kms,jm2)
1729	CALL wrf_debug(200, message)
1730	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1731	CALL wrf_debug(200, message)
1732	WRITE(message,*)' al: ',grid%al(im2,km2,jm2), ' p: ', grid%p(im2,1,jm2), &
1733	' ph: ',grid%ph_2(im2,1,jm2)
1734	CALL wrf_debug(200, message)
1735	WRITE(message,*)' mu: ',grid%mu_2(im2,jm2), 'alt: ', grid%alt(im2,1,jm2), &
1736	' mu: ',grid%mu_2(im2,jm2), ' znu: ', grid%znu(1),' ph 1: ',grid%ph_2(im2,2,jm2)
1737	CALL wrf_debug(200, message)
1738	PRINT *,' c2a: ',c2a(im2,1,jm2), 't_2: ', grid%t_2(im2,1,jm2), &
1739	' t_save: ',grid%t_save(im2,1,jm2),' pm1: ',pm1(im2,1,jm2)
1740	CALL wrf_debug(200, message)
1741	#endif
1742
1743	!-----------------------------------------------------------
1744	! acoustic integration polar filter for smallstep w, geopotential
1745	!-----------------------------------------------------------
1746
1747	IF ( (config_flags%polar) .AND. (config_flags%non_hydrostatic) ) THEN
1748
1749	CALL pxft ( grid=grid &
1750	,lineno=__LINE__ &
1751	,flag_uv = 0 &
1752	,flag_rurv = 0 &
1753	,flag_wph = 1 &
1754	,flag_ww = 0 &
1755	,flag_t = 0 &
1756	,flag_mu = 0 &
1757	,flag_mut = 0 &
1758	,flag_moist = 0 &
1759	,flag_chem = 0 &
1760	,flag_tracer = 0 &
1761	,flag_scalar = 0 &
1762	,positive_definite = .FALSE. &
1763	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
1764	,fft_filter_lat = config_flags%fft_filter_lat &
1765	,dclat = dclat &
1766	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
1767	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
1768	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
1769	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
1770	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
1771
1772	END IF
1773
1774	!-----------------------------------------------------------
1775	! end acoustic integration polar filter for smallstep w, geopotential
1776	!-----------------------------------------------------------
1777
1778	!$OMP PARALLEL DO &
1779	!$OMP PRIVATE ( ij )
1780	DO ij = 1 , grid%num_tiles
1781
1782	BENCH_START(sumflux_tim)
1783	CALL sumflux ( grid%u_2, grid%v_2, grid%ww, &
1784	grid%u_save, grid%v_save, ww1, &
1785	grid%muu, grid%muv, &
1786	grid%ru_m, grid%rv_m, grid%ww_m, grid%epssm, &
1787	grid%msfux, grid% msfuy, grid%msfvx, &
1788	grid%msfvx_inv, grid%msfvy, &
1789	iteration, number_of_small_timesteps, &
1790	ids, ide, jds, jde, kds, kde, &
1791	ims, ime, jms, jme, kms, kme, &
1792	grid%i_start(ij), grid%i_end(ij), &
1793	grid%j_start(ij), grid%j_end(ij), &
1794	k_start , k_end )
1795	BENCH_END(sumflux_tim)
1796	#ifdef LMDZ1
1797	WRITE(message, *)' dyn_em: after sumflux'
1798	CALL wrf_debug(200, message)
1799	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1800	' u_tend: ', ru_tendf(im2,1,jm2)
1801	CALL wrf_debug(200, message)
1802	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1803	'p sfc: ',p8w(im2,kms,jm2)
1804	CALL wrf_debug(200, message)
1805	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1806	CALL wrf_debug(200, message)
1807	WRITE(message,*)' al: ',grid%al(im2,km2,jm2), ' p: ', grid%p(im2,1,jm2), &
1808	' ph: ',grid%ph_2(im2,1,jm2)
1809	CALL wrf_debug(200, message)
1810	WRITE(message,*)' mu: ',grid%mu_2(im2,jm2), 'alt: ', grid%alt(im2,1,jm2), &
1811	' mu: ',grid%mu_2(im2,jm2), ' znu: ', grid%znu(1),' ph 1: ',grid%ph_2(im2,2,jm2)
1812	CALL wrf_debug(200, message)
1813	PRINT *,' c2a: ',c2a(im2,1,jm2), 't_2: ', grid%t_2(im2,1,jm2), &
1814	' t_save: ',grid%t_save(im2,1,jm2),' pm1: ',pm1(im2,1,jm2)
1815	CALL wrf_debug(200, message)
1816	#endif
1817
1818	IF( config_flags%specified .or. config_flags%nested ) THEN
1819
1820	BENCH_START(spec_bdynhyd_tim)
1821	IF (config_flags%non_hydrostatic) THEN
1822	CALL spec_bdyupdate_ph( ph_save, grid%ph_2, ph_tend, &
1823	mu_tend, grid%muts, dts_rk, &
1824	'h' , config_flags, &
1825	grid%spec_zone, &
1826	ids,ide, jds,jde, kds,kde, &
1827	ims,ime, jms,jme, kms,kme, &
1828	ips,ipe, jps,jpe, kps,kpe, &
1829	grid%i_start(ij), grid%i_end(ij),&
1830	grid%j_start(ij), grid%j_end(ij),&
1831	k_start , k_end )
1832	#ifdef LMDZ1
1833	WRITE(message, *)' dyn_em: after spec_bdynhyd_ph'
1834	CALL wrf_debug(200, message)
1835	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1836	' u_tend: ', ru_tendf(im2,1,jm2)
1837	CALL wrf_debug(200, message)
1838	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1839	'p sfc: ',p8w(im2,kms,jm2)
1840	CALL wrf_debug(200, message)
1841	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1842	CALL wrf_debug(200, message)
1843	WRITE(message,*)' al: ',grid%al(im2,km2,jm2), ' p: ', grid%p(im2,1,jm2), &
1844	' ph: ',grid%ph_2(im2,1,jm2)
1845	CALL wrf_debug(200, message)
1846	WRITE(message,*)' mu: ',grid%mu_2(im2,jm2), 'alt: ', grid%alt(im2,1,jm2), &
1847	' mu: ',grid%mu_2(im2,jm2), ' znu: ', grid%znu(1),' ph 1: ',grid%ph_2(im2,2,jm2)
1848	CALL wrf_debug(200, message)
1849	PRINT *,' c2a: ',c2a(im2,1,jm2), 't_2: ', grid%t_2(im2,1,jm2), &
1850	' t_save: ',grid%t_save(im2,1,jm2),' pm1: ',pm1(im2,1,jm2)
1851	CALL wrf_debug(200, message)
1852	#endif
1853	IF( config_flags%specified ) THEN
1854	CALL zero_grad_bdy ( grid%w_2, &
1855	'w' , config_flags, &
1856	grid%spec_zone, &
1857	ids,ide, jds,jde, kds,kde, &
1858	ims,ime, jms,jme, kms,kme, &
1859	ips,ipe, jps,jpe, kps,kpe, &
1860	grid%i_start(ij), grid%i_end(ij), &
1861	grid%j_start(ij), grid%j_end(ij), &
1862	k_start , k_end )
1863	ELSE
1864	CALL spec_bdyupdate ( grid%w_2, rw_tend, dts_rk, &
1865	'h' , config_flags, &
1866	grid%spec_zone, &
1867	ids,ide, jds,jde, kds,kde, &
1868	ims,ime, jms,jme, kms,kme, &
1869	ips,ipe, jps,jpe, kps,kpe, &
1870	grid%i_start(ij), grid%i_end(ij),&
1871	grid%j_start(ij), grid%j_end(ij),&
1872	k_start , k_end )
1873	ENDIF
1874	ENDIF
1875	BENCH_END(spec_bdynhyd_tim)
1876	ENDIF
1877	#ifdef LMDZ1
1878	WRITE(message, *)' dyn_em: after spec_bdynhyd_tim'
1879	CALL wrf_debug(200, message)
1880	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1881	' u_tend: ', ru_tendf(im2,1,jm2)
1882	CALL wrf_debug(200, message)
1883	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1884	'p sfc: ',p8w(im2,kms,jm2)
1885	CALL wrf_debug(200, message)
1886	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1887	CALL wrf_debug(200, message)
1888	WRITE(message,*)' al: ',grid%al(im2,km2,jm2), ' p: ', grid%p(im2,1,jm2), &
1889	' ph: ',grid%ph_2(im2,1,jm2)
1890	CALL wrf_debug(200, message)
1891	WRITE(message,*)' mu: ',grid%mu_2(im2,jm2), 'alt: ', grid%alt(im2,1,jm2), &
1892	' mu: ',grid%mu_2(im2,jm2), ' znu: ', grid%znu(1),' ph 1: ',grid%ph_2(im2,2,jm2)
1893	CALL wrf_debug(200, message)
1894	PRINT *,' c2a: ',c2a(im2,1,jm2), 't_2: ', grid%t_2(im2,1,jm2), &
1895	' t_save: ',grid%t_save(im2,1,jm2),' pm1: ',pm1(im2,1,jm2)
1896	CALL wrf_debug(200, message)
1897	#endif
1898
1899	BENCH_START(cald_p_rho_tim)
1900	CALL calc_p_rho( grid%al, grid%p, grid%ph_2, &
1901	grid%alt, grid%t_2, grid%t_save, c2a, pm1, &
1902	grid%mu_2, grid%muts, grid%znu, t0, &
1903	grid%rdnw, grid%dnw, grid%smdiv, &
1904	config_flags%non_hydrostatic, iteration, &
1905	ids, ide, jds, jde, kds, kde, &
1906	ims, ime, jms, jme, kms, kme, &
1907	grid%i_start(ij), grid%i_end(ij), &
1908	grid%j_start(ij), grid%j_end(ij), &
1909	k_start , k_end )
1910	BENCH_END(cald_p_rho_tim)
1911	#ifdef LMDZ1
1912	WRITE(message, *)' dyn_em: after calc_p_rho'
1913	CALL wrf_debug(200, message)
1914	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1915	' u_tend: ', ru_tendf(im2,1,jm2)
1916	CALL wrf_debug(200, message)
1917	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1918	'p sfc: ',p8w(im2,kms,jm2)
1919	CALL wrf_debug(200, message)
1920	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1921	CALL wrf_debug(200, message)
1922	#endif
1923
1924	ENDDO
1925	!$OMP END PARALLEL DO
1926	#ifdef LMDZ1
1927	WRITE(message, *)' dyn_em: after geopotential'
1928	CALL wrf_debug(200, message)
1929	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1930	' u_tend: ', ru_tendf(im2,1,jm2)
1931	CALL wrf_debug(200, message)
1932	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1933	'p sfc: ',p8w(im2,kms,jm2)
1934	CALL wrf_debug(200, message)
1935	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1936	CALL wrf_debug(200, message)
1937	#endif
1938
1939	#ifdef DM_PARALLEL
1940	!
1941	! Stencils for patch communications (WCS, 29 June 2001)
1942	!
1943	! * *
1944	! * + * * + * +
1945	! * *
1946	!
1947	! ph_2 x
1948	! al x
1949	! p x
1950	!
1951	! 2D variables (x,y)
1952	!
1953	! mu_2 x
1954	! muts x
1955	! mudf x
1956
1957	# include "HALO_EM_C2.inc"
1958	# include "PERIOD_BDY_EM_B3.inc"
1959	#endif
1960
1961	BENCH_START(phys_bc_tim)
1962	!$OMP PARALLEL DO &
1963	!$OMP PRIVATE ( ij )
1964	DO ij = 1 , grid%num_tiles
1965
1966	! boundary condition set for next small timestep
1967
1968	CALL set_physical_bc3d( grid%ph_2, 'w', config_flags, &
1969	ids, ide, jds, jde, kds, kde, &
1970	ims, ime, jms, jme, kms, kme, &
1971	ips, ipe, jps, jpe, kps, kpe, &
1972	grid%i_start(ij), grid%i_end(ij), &
1973	grid%j_start(ij), grid%j_end(ij), &
1974	k_start , k_end )
1975
1976	CALL set_physical_bc3d( grid%al, 'p', config_flags, &
1977	ids, ide, jds, jde, kds, kde, &
1978	ims, ime, jms, jme, kms, kme, &
1979	ips, ipe, jps, jpe, kps, kpe, &
1980	grid%i_start(ij), grid%i_end(ij), &
1981	grid%j_start(ij), grid%j_end(ij), &
1982	k_start , k_end )
1983
1984	CALL set_physical_bc3d( grid%p, 'p', config_flags, &
1985	ids, ide, jds, jde, kds, kde, &
1986	ims, ime, jms, jme, kms, kme, &
1987	ips, ipe, jps, jpe, kps, kpe, &
1988	grid%i_start(ij), grid%i_end(ij), &
1989	grid%j_start(ij), grid%j_end(ij), &
1990	k_start , k_end )
1991
1992	CALL set_physical_bc2d( grid%muts, 't', config_flags, &
1993	ids, ide, jds, jde, &
1994	ims, ime, jms, jme, &
1995	ips, ipe, jps, jpe, &
1996	grid%i_start(ij), grid%i_end(ij), &
1997	grid%j_start(ij), grid%j_end(ij) )
1998
1999	CALL set_physical_bc2d( grid%mu_2, 't', config_flags, &
2000	ids, ide, jds, jde, &
2001	ims, ime, jms, jme, &
2002	ips, ipe, jps, jpe, &
2003	grid%i_start(ij), grid%i_end(ij), &
2004	grid%j_start(ij), grid%j_end(ij) )
2005
2006	CALL set_physical_bc2d( grid%mudf, 't', config_flags, &
2007	ids, ide, jds, jde, &
2008	ims, ime, jms, jme, &
2009	ips, ipe, jps, jpe, &
2010	grid%i_start(ij), grid%i_end(ij), &
2011	grid%j_start(ij), grid%j_end(ij) )
2012
2013	END DO
2014	!$OMP END PARALLEL DO
2015	BENCH_END(phys_bc_tim)
2016
2017	END DO small_steps
2018
2019	!$OMP PARALLEL DO &
2020	!$OMP PRIVATE ( ij )
2021	DO ij = 1 , grid%num_tiles
2022
2023	CALL wrf_debug ( 200 , ' call rk_small_finish' )
2024
2025	! change time-perturbation variables back to
2026	! full perturbation variables.
2027	! first get updated mu at u and v points
2028
2029	BENCH_START(calc_mu_uv_tim)
2030	CALL calc_mu_uv_1 ( config_flags, &
2031	grid%muts, muus, muvs, &
2032	ids, ide, jds, jde, kds, kde, &
2033	ims, ime, jms, jme, kms, kme, &
2034	grid%i_start(ij), grid%i_end(ij), &
2035	grid%j_start(ij), grid%j_end(ij), &
2036	k_start , k_end )
2037	BENCH_END(calc_mu_uv_tim)
2038	BENCH_START(small_step_finish_tim)
2039	CALL small_step_finish( grid%u_2, grid%u_1, grid%v_2, grid%v_1, grid%w_2, grid%w_1, &
2040	grid%t_2, grid%t_1, grid%ph_2, grid%ph_1, grid%ww, ww1, &
2041	grid%mu_2, grid%mu_1, &
2042	grid%mut, grid%muts, grid%muu, muus, grid%muv, muvs, &
2043	grid%u_save, grid%v_save, w_save, &
2044	grid%t_save, ph_save, mu_save, &
2045	grid%msfux,grid%msfuy, grid%msfvx,grid%msfvy, grid%msftx,grid%msfty, &
2046	grid%h_diabatic, &
2047	number_of_small_timesteps,dts_rk, &
2048	rk_step, rk_order, &
2049	ids, ide, jds, jde, kds, kde, &
2050	ims, ime, jms, jme, kms, kme, &
2051	grid%i_start(ij), grid%i_end(ij), &
2052	grid%j_start(ij), grid%j_end(ij), &
2053	k_start , k_end )
2054	! call to set ru_m, rv_m and ww_m b.c's for PD advection
2055
2056	IF (rk_step == rk_order) THEN
2057
2058	CALL set_physical_bc3d( grid%ru_m, 'u', config_flags, &
2059	ids, ide, jds, jde, kds, kde, &
2060	ims, ime, jms, jme, kms, kme, &
2061	ips, ipe, jps, jpe, kps, kpe, &
2062	grid%i_start(ij), grid%i_end(ij), &
2063	grid%j_start(ij), grid%j_end(ij), &
2064	k_start , k_end )
2065
2066	CALL set_physical_bc3d( grid%rv_m, 'v', config_flags, &
2067	ids, ide, jds, jde, kds, kde, &
2068	ims, ime, jms, jme, kms, kme, &
2069	ips, ipe, jps, jpe, kps, kpe, &
2070	grid%i_start(ij), grid%i_end(ij), &
2071	grid%j_start(ij), grid%j_end(ij), &
2072	k_start , k_end )
2073
2074	CALL set_physical_bc3d( grid%ww_m, 'w', config_flags, &
2075	ids, ide, jds, jde, kds, kde, &
2076	ims, ime, jms, jme, kms, kme, &
2077	ips, ipe, jps, jpe, kps, kpe, &
2078	grid%i_start(ij), grid%i_end(ij), &
2079	grid%j_start(ij), grid%j_end(ij), &
2080	k_start , k_end )
2081
2082	CALL set_physical_bc2d( grid%mut, 't', config_flags, &
2083	ids, ide, jds, jde, &
2084	ims, ime, jms, jme, &
2085	ips, ipe, jps, jpe, &
2086	grid%i_start(ij), grid%i_end(ij), &
2087	grid%j_start(ij), grid%j_end(ij) )
2088
2089	CALL set_physical_bc2d( grid%muts, 't', config_flags, &
2090	ids, ide, jds, jde, &
2091	ims, ime, jms, jme, &
2092	ips, ipe, jps, jpe, &
2093	grid%i_start(ij), grid%i_end(ij), &
2094	grid%j_start(ij), grid%j_end(ij) )
2095
2096	END IF
2097
2098	BENCH_END(small_step_finish_tim)
2099
2100	END DO
2101	!$OMP END PARALLEL DO
2102	#ifdef LMDZ1
2103	WRITE(message, *)' dyn_em: after rk_small_finish'
2104	CALL wrf_debug(200, message)
2105	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
2106	' u_tend: ', ru_tendf(im2,1,jm2)
2107	CALL wrf_debug(200, message)
2108	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
2109	'p sfc: ',p8w(im2,kms,jm2)
2110	CALL wrf_debug(200, message)
2111	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
2112	CALL wrf_debug(200, message)
2113	#endif
2114
2115	!-----------------------------------------------------------
2116	! polar filter for full dynamics variables and time-averaged mass fluxes
2117	!-----------------------------------------------------------
2118
2119	IF (config_flags%polar) THEN
2120
2121	CALL pxft ( grid=grid &
2122	,lineno=__LINE__ &
2123	,flag_uv = 1 &
2124	,flag_rurv = 1 &
2125	,flag_wph = 1 &
2126	,flag_ww = 1 &
2127	,flag_t = 1 &
2128	,flag_mu = 1 &
2129	,flag_mut = 1 &
2130	,flag_moist = 0 &
2131	,flag_chem = 0 &
2132	,flag_tracer = 0 &
2133	,flag_scalar = 0 &
2134	,positive_definite = .FALSE. &
2135	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
2136	,fft_filter_lat = config_flags%fft_filter_lat &
2137	,dclat = dclat &
2138	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
2139	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
2140	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
2141	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
2142	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
2143
2144	END IF
2145
2146	!-----------------------------------------------------------
2147	! end polar filter for full dynamics variables and time-averaged mass fluxes
2148	!-----------------------------------------------------------
2149
2150	!-----------------------------------------------------------------------
2151	! add in physics tendency first if positive definite advection is used.
2152	! pd advection applies advective flux limiter on last runge-kutta step
2153	!-----------------------------------------------------------------------
2154	! first moisture
2155
2156	IF ((config_flags%moist_adv_opt /= ORIGINAL) .and. (rk_step == rk_order)) THEN
2157
2158	!$OMP PARALLEL DO &
2159	!$OMP PRIVATE ( ij )
2160	DO ij = 1 , grid%num_tiles
2161	CALL wrf_debug ( 200 , ' call rk_update_scalar_pd' )
2162	DO im = PARAM_FIRST_SCALAR, num_3d_m
2163	CALL rk_update_scalar_pd( im, im, &
2164	moist_old(ims,kms,jms,im), &
2165	moist_tend(ims,kms,jms,im), &
2166	grid%mu_1, grid%mu_1, grid%mub, &
2167	rk_step, dt_rk, grid%spec_zone, &
2168	config_flags, &
2169	ids, ide, jds, jde, kds, kde, &
2170	ims, ime, jms, jme, kms, kme, &
2171	grid%i_start(ij), grid%i_end(ij), &
2172	grid%j_start(ij), grid%j_end(ij), &
2173	k_start , k_end )
2174	ENDDO
2175	END DO
2176	!$OMP END PARALLEL DO
2177
2178	!---------------------- positive definite bc call
2179	#ifdef DM_PARALLEL
2180	IF (config_flags%moist_adv_opt /= ORIGINAL) THEN
2181	IF ( config_flags%h_sca_adv_order <= 4 ) THEN
2182	# include "HALO_EM_MOIST_OLD_E_5.inc"
2183	ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN
2184	# include "HALO_EM_MOIST_OLD_E_7.inc"
2185	ELSE
2186	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
2187	CALL wrf_error_fatal(TRIM(wrf_err_message))
2188	ENDIF
2189	ENDIF
2190	#endif
2191
2192	#ifdef DM_PARALLEL
2193	# include "PERIOD_BDY_EM_MOIST_OLD.inc"
2194	#endif
2195
2196	!$OMP PARALLEL DO &
2197	!$OMP PRIVATE ( ij )
2198	DO ij = 1 , grid%num_tiles
2199	IF (num_3d_m >= PARAM_FIRST_SCALAR) THEN
2200	DO im = PARAM_FIRST_SCALAR , num_3d_m
2201	CALL set_physical_bc3d( moist_old(ims,kms,jms,im), 'p', config_flags, &
2202	ids, ide, jds, jde, kds, kde, &
2203	ims, ime, jms, jme, kms, kme, &
2204	ips, ipe, jps, jpe, kps, kpe, &
2205	grid%i_start(ij), grid%i_end(ij), &
2206	grid%j_start(ij), grid%j_end(ij), &
2207	k_start , k_end )
2208	END DO
2209	ENDIF
2210	END DO
2211	!$OMP END PARALLEL DO
2212
2213	END IF ! end if for moist_adv_opt
2214
2215	! scalars
2216
2217	IF ((config_flags%scalar_adv_opt /= ORIGINAL) .and. (rk_step == rk_order)) THEN
2218
2219	!$OMP PARALLEL DO &
2220	!$OMP PRIVATE ( ij )
2221	DO ij = 1 , grid%num_tiles
2222	CALL wrf_debug ( 200 , ' call rk_update_scalar_pd' )
2223	DO im = PARAM_FIRST_SCALAR, num_3d_s
2224	CALL rk_update_scalar_pd( im, im, &
2225	scalar_old(ims,kms,jms,im), &
2226	scalar_tend(ims,kms,jms,im), &
2227	grid%mu_1, grid%mu_1, grid%mub, &
2228	rk_step, dt_rk, grid%spec_zone, &
2229	config_flags, &
2230	ids, ide, jds, jde, kds, kde, &
2231	ims, ime, jms, jme, kms, kme, &
2232	grid%i_start(ij), grid%i_end(ij), &
2233	grid%j_start(ij), grid%j_end(ij), &
2234	k_start , k_end )
2235	ENDDO
2236	ENDDO
2237	!$OMP END PARALLEL DO
2238
2239	!---------------------- positive definite bc call
2240	#ifdef DM_PARALLEL
2241	IF (config_flags%scalar_adv_opt /= ORIGINAL) THEN
2242	#ifndef RSL
2243	IF ( config_flags%h_sca_adv_order <= 4 ) THEN
2244	# include "HALO_EM_SCALAR_OLD_E_5.inc"
2245	ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN
2246	# include "HALO_EM_SCALAR_OLD_E_7.inc"
2247	ELSE
2248	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
2249	CALL wrf_error_fatal(TRIM(wrf_err_message))
2250	ENDIF
2251	#else
2252	WRITE(wrf_err_message,*)'cannot use pd scheme with RSL - use RSL-LITE'
2253	CALL wrf_error_fatal(TRIM(wrf_err_message))
2254	#endif
2255	endif
2256	#endif
2257
2258	#ifdef DM_PARALLEL
2259	# include "PERIOD_BDY_EM_SCALAR_OLD.inc"
2260	#endif
2261	!$OMP PARALLEL DO &
2262	!$OMP PRIVATE ( ij )
2263
2264	DO ij = 1 , grid%num_tiles
2265	IF (num_3d_s >= PARAM_FIRST_SCALAR) THEN
2266	DO im = PARAM_FIRST_SCALAR , num_3d_s
2267	CALL set_physical_bc3d( scalar_old(ims,kms,jms,im), 'p', config_flags, &
2268	ids, ide, jds, jde, kds, kde, &
2269	ims, ime, jms, jme, kms, kme, &
2270	ips, ipe, jps, jpe, kps, kpe, &
2271	grid%i_start(ij), grid%i_end(ij), &
2272	grid%j_start(ij), grid%j_end(ij), &
2273	k_start , k_end )
2274	END DO
2275	ENDIF
2276	END DO
2277	!$OMP END PARALLEL DO
2278
2279	END IF ! end if for scalar_adv_opt
2280
2281	! chem
2282
2283	IF ((config_flags%chem_adv_opt /= ORIGINAL) .and. (rk_step == rk_order)) THEN
2284
2285	!$OMP PARALLEL DO &
2286	!$OMP PRIVATE ( ij )
2287	DO ij = 1 , grid%num_tiles
2288	CALL wrf_debug ( 200 , ' call rk_update_scalar_pd' )
2289	DO im = PARAM_FIRST_SCALAR, num_3d_c
2290	CALL rk_update_scalar_pd( im, im, &
2291	chem_old(ims,kms,jms,im), &
2292	chem_tend(ims,kms,jms,im), &
2293	grid%mu_1, grid%mu_1, grid%mub, &
2294	rk_step, dt_rk, grid%spec_zone, &
2295	config_flags, &
2296	ids, ide, jds, jde, kds, kde, &
2297	ims, ime, jms, jme, kms, kme, &
2298	grid%i_start(ij), grid%i_end(ij), &
2299	grid%j_start(ij), grid%j_end(ij), &
2300	k_start , k_end )
2301	ENDDO
2302	END DO
2303	!$OMP END PARALLEL DO
2304
2305	!---------------------- positive definite bc call
2306	#ifdef DM_PARALLEL
2307	IF (config_flags%chem_adv_opt /= ORIGINAL) THEN
2308	IF ( config_flags%h_sca_adv_order <= 4 ) THEN
2309	# include "HALO_EM_CHEM_OLD_E_5.inc"
2310	ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN
2311	# include "HALO_EM_CHEM_OLD_E_7.inc"
2312	ELSE
2313	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
2314	CALL wrf_error_fatal(TRIM(wrf_err_message))
2315	ENDIF
2316	ENDIF
2317	#endif
2318
2319	#ifdef DM_PARALLEL
2320	# include "PERIOD_BDY_EM_CHEM_OLD.inc"
2321	#endif
2322
2323	!$OMP PARALLEL DO &
2324	!$OMP PRIVATE ( ij )
2325	DO ij = 1 , grid%num_tiles
2326	IF (num_3d_c >= PARAM_FIRST_SCALAR) THEN
2327	DO im = PARAM_FIRST_SCALAR , num_3d_c
2328	CALL set_physical_bc3d( chem_old(ims,kms,jms,im), 'p', config_flags, &
2329	ids, ide, jds, jde, kds, kde, &
2330	ims, ime, jms, jme, kms, kme, &
2331	ips, ipe, jps, jpe, kps, kpe, &
2332	grid%i_start(ij), grid%i_end(ij), &
2333	grid%j_start(ij), grid%j_end(ij), &
2334	k_start , k_end )
2335	END DO
2336	ENDIF
2337	END DO
2338	!$OMP END PARALLEL DO
2339
2340	ENDIF ! end if for chem_adv_opt
2341
2342	! tracer
2343
2344	IF ((config_flags%tracer_adv_opt /= ORIGINAL) .and. (rk_step == rk_order)) THEN
2345
2346	!$OMP PARALLEL DO &
2347	!$OMP PRIVATE ( ij )
2348	DO ij = 1 , grid%num_tiles
2349	CALL wrf_debug ( 200 , ' call rk_update_scalar_pd' )
2350	DO im = PARAM_FIRST_SCALAR, num_tracer
2351	CALL rk_update_scalar_pd( im, im, &
2352	tracer_old(ims,kms,jms,im), &
2353	tracer_tend(ims,kms,jms,im), &
2354	grid%mu_1, grid%mu_1, grid%mub, &
2355	rk_step, dt_rk, grid%spec_zone, &
2356	config_flags, &
2357	ids, ide, jds, jde, kds, kde, &
2358	ims, ime, jms, jme, kms, kme, &
2359	grid%i_start(ij), grid%i_end(ij), &
2360	grid%j_start(ij), grid%j_end(ij), &
2361	k_start , k_end )
2362	ENDDO
2363	END DO
2364	!$OMP END PARALLEL DO
2365
2366	!---------------------- positive definite bc call
2367	#ifdef DM_PARALLEL
2368	IF (config_flags%tracer_adv_opt /= ORIGINAL) THEN
2369	IF ( config_flags%h_sca_adv_order <= 4 ) THEN
2370	# include "HALO_EM_TRACER_OLD_E_5.inc"
2371	ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN
2372	# include "HALO_EM_TRACER_OLD_E_7.inc"
2373	ELSE
2374	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
2375	CALL wrf_error_fatal(TRIM(wrf_err_message))
2376	ENDIF
2377	ENDIF
2378	#endif
2379
2380	#ifdef DM_PARALLEL
2381	# include "PERIOD_BDY_EM_TRACER_OLD.inc"
2382	#endif
2383
2384	!$OMP PARALLEL DO &
2385	!$OMP PRIVATE ( ij )
2386	DO ij = 1 , grid%num_tiles
2387	IF (num_tracer >= PARAM_FIRST_SCALAR) THEN
2388	DO im = PARAM_FIRST_SCALAR , num_tracer
2389	CALL set_physical_bc3d( tracer_old(ims,kms,jms,im), 'p', config_flags, &
2390	ids, ide, jds, jde, kds, kde, &
2391	ims, ime, jms, jme, kms, kme, &
2392	ips, ipe, jps, jpe, kps, kpe, &
2393	grid%i_start(ij), grid%i_end(ij), &
2394	grid%j_start(ij), grid%j_end(ij), &
2395	k_start , k_end )
2396	END DO
2397	ENDIF
2398	END DO
2399	!$OMP END PARALLEL DO
2400
2401	ENDIF ! end if for tracer_adv_opt
2402
2403	! tke
2404
2405	IF ((config_flags%tke_adv_opt /= ORIGINAL) .and. (rk_step == rk_order) &
2406	.and. (config_flags%km_opt .eq. 2) ) THEN
2407
2408	!$OMP PARALLEL DO &
2409	!$OMP PRIVATE ( ij )
2410	DO ij = 1 , grid%num_tiles
2411	CALL wrf_debug ( 200 , ' call rk_update_scalar_pd' )
2412	CALL rk_update_scalar_pd( 1, 1, &
2413	grid%tke_1, &
2414	tke_tend(ims,kms,jms), &
2415	grid%mu_1, grid%mu_1, grid%mub, &
2416	rk_step, dt_rk, grid%spec_zone, &
2417	config_flags, &
2418	ids, ide, jds, jde, kds, kde, &
2419	ims, ime, jms, jme, kms, kme, &
2420	grid%i_start(ij), grid%i_end(ij), &
2421	grid%j_start(ij), grid%j_end(ij), &
2422	k_start , k_end )
2423	ENDDO
2424	!$OMP END PARALLEL DO
2425
2426	!---------------------- positive definite bc call
2427	#ifdef DM_PARALLEL
2428	IF (config_flags%tke_adv_opt /= ORIGINAL) THEN
2429	IF ( config_flags%h_sca_adv_order <= 4 ) THEN
2430	# include "HALO_EM_TKE_OLD_E_5.inc"
2431	ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN
2432	# include "HALO_EM_TKE_OLD_E_7.inc"
2433	ELSE
2434	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
2435	CALL wrf_error_fatal(TRIM(wrf_err_message))
2436	ENDIF
2437	ENDIF
2438	#endif
2439
2440	#ifdef DM_PARALLEL
2441	# include "PERIOD_BDY_EM_TKE_OLD.inc"
2442	#endif
2443
2444	!$OMP PARALLEL DO &
2445	!$OMP PRIVATE ( ij )
2446	DO ij = 1 , grid%num_tiles
2447	CALL set_physical_bc3d( grid%tke_1, 'p', config_flags, &
2448	ids, ide, jds, jde, kds, kde, &
2449	ims, ime, jms, jme, kms, kme, &
2450	ips, ipe, jps, jpe, kps, kpe, &
2451	grid%i_start(ij), grid%i_end(ij), &
2452	grid%j_start(ij), grid%j_end(ij), &
2453	k_start , k_end )
2454	END DO
2455	!$OMP END PARALLEL DO
2456
2457	!--- end of positive definite physics tendency update
2458
2459	END IF ! end if for tke_adv_opt
2460
2461	#ifdef DM_PARALLEL
2462	!
2463	! Stencils for patch communications (WCS, 29 June 2001)
2464	!
2465	! * * * * *
2466	! * * * * *
2467	! * * + * *
2468	! * * * * *
2469	! * * * * *
2470	!
2471	! ru_m x
2472	! rv_m x
2473	! ww_m x
2474	! mut x
2475	!
2476	!--------------------------------------------------------------
2477
2478	# include "HALO_EM_D.inc"
2479	! WCS addition 11/19/08
2480	# include "PERIOD_EM_DA.inc"
2481	#endif
2482
2483	!<DESCRIPTION>
2484	!<pre>
2485	! (4) Still within the RK loop, the scalar variables are advanced.
2486	!
2487	! For the moist and chem variables, each one is advanced
2488	! individually, using named loops "moist_variable_loop:"
2489	! and "chem_variable_loop:". Each RK substep begins by
2490	! calculating the advective tendency, and, for the first RK step,
2491	! 3D mixing (calling rk_scalar_tend) followed by an update
2492	! of the scalar (calling rk_update_scalar).
2493	!</pre>
2494	!</DESCRIPTION>
2495
2496	PRINT *,' Lluis PARAM_FIRST_SCALAR: ', PARAM_FIRST_SCALAR, &
2497	' p_qv: ',p_qv, ' p_qc: ',p_qc,' U(moist): ', UBOUND(moist)
2498	PRINT *,' moist: ',moist(config_flags%i_check_point, 2, &
2499	config_flags%j_check_point,:)
2500
2501	moist_scalar_advance: IF (num_3d_m >= PARAM_FIRST_SCALAR ) THEN
2502
2503	moist_variable_loop: DO im = PARAM_FIRST_SCALAR, num_3d_m
2504
2505	! adv_moist_cond is set in module_physics_init based on mp_physics choice
2506	! true except for Ferrier scheme
2507
2508	IF (grid%adv_moist_cond .or. im==p_qv ) THEN
2509
2510	!$OMP PARALLEL DO &
2511	!$OMP PRIVATE ( ij )
2512	moist_tile_loop_1: DO ij = 1 , grid%num_tiles
2513
2514	CALL wrf_debug ( 200 , ' call rk_scalar_tend' )
2515	tenddec = .false.
2516	PRINT *,' Lluis before rk_update_scalar moist im=',im, &
2517	' moist_old moist_tend moist________'
2518	DO iz = kms, kme
2519	PRINT *,moist_old(config_flags%i_check_point, iz, &
2520	config_flags%j_check_point,im), &
2521	moist_tend(config_flags%i_check_point, iz, &
2522	config_flags%j_check_point,im), &
2523	moist(config_flags%i_check_point, iz, &
2524	config_flags%j_check_point,im)
2525	END DO
2526
2527	BENCH_START(rk_scalar_tend_tim)
2528	CALL rk_scalar_tend ( im, im, config_flags, tenddec, &
2529	rk_step, dt_rk, &
2530	grid%ru_m, grid%rv_m, grid%ww_m, &
2531	grid%muts, grid%mub, grid%mu_1, &
2532	grid%alt, &
2533	moist_old(ims,kms,jms,im), &
2534	moist(ims,kms,jms,im), &
2535	moist_tend(ims,kms,jms,im), &
2536	advect_tend,h_tendency,z_tendency,grid%rqvften, &
2537	grid%qv_base, .true., grid%fnm, grid%fnp, &
2538	grid%msfux,grid%msfuy, grid%msfvx, grid%msfvx_inv,&
2539	grid%msfvy, grid%msftx,grid%msfty, &
2540	grid%rdx, grid%rdy, grid%rdn, grid%rdnw, grid%khdif, &
2541	grid%kvdif, grid%xkhh, &
2542	grid%diff_6th_opt, grid%diff_6th_factor, &
2543	config_flags%moist_adv_opt, &
2544	ids, ide, jds, jde, kds, kde, &
2545	ims, ime, jms, jme, kms, kme, &
2546	grid%i_start(ij), grid%i_end(ij), &
2547	grid%j_start(ij), grid%j_end(ij), &
2548	k_start , k_end )
2549
2550	BENCH_END(rk_scalar_tend_tim)
2551
2552	BENCH_START(rlx_bdy_scalar_tim)
2553	IF( ( config_flags%specified .or. config_flags%nested ) .and. rk_step == 1 ) THEN
2554	IF ( im .EQ. P_QV .OR. config_flags%nested ) THEN
2555	CALL relax_bdy_scalar ( moist_tend(ims,kms,jms,im), &
2556	moist(ims,kms,jms,im), grid%mut, &
2557	moist_bxs(jms,kms,1,im),moist_bxe(jms,kms,1,im), &
2558	moist_bys(ims,kms,1,im),moist_bye(ims,kms,1,im), &
2559	moist_btxs(jms,kms,1,im),moist_btxe(jms,kms,1,im), &
2560	moist_btys(ims,kms,1,im),moist_btye(ims,kms,1,im), &
2561	config_flags%spec_bdy_width, grid%spec_zone, grid%relax_zone, &
2562	grid%dtbc, grid%fcx, grid%gcx, &
2563	config_flags, &
2564	ids,ide, jds,jde, kds,kde, & ! domain dims
2565	ims,ime, jms,jme, kms,kme, & ! memory dims
2566	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
2567	grid%i_start(ij), grid%i_end(ij), &
2568	grid%j_start(ij), grid%j_end(ij), &
2569	k_start, k_end )
2570
2571	CALL spec_bdy_scalar ( moist_tend(ims,kms,jms,im), &
2572	moist_bxs(jms,kms,1,im),moist_bxe(jms,kms,1,im), &
2573	moist_bys(ims,kms,1,im),moist_bye(ims,kms,1,im), &
2574	moist_btxs(jms,kms,1,im),moist_btxe(jms,kms,1,im), &
2575	moist_btys(ims,kms,1,im),moist_btye(ims,kms,1,im), &
2576	config_flags%spec_bdy_width, grid%spec_zone, &
2577	config_flags, &
2578	ids,ide, jds,jde, kds,kde, & ! domain dims
2579	ims,ime, jms,jme, kms,kme, & ! memory dims
2580	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
2581	grid%i_start(ij), grid%i_end(ij), &
2582	grid%j_start(ij), grid%j_end(ij), &
2583	k_start, k_end )
2584	ENDIF
2585	ENDIF
2586	BENCH_END(rlx_bdy_scalar_tim)
2587
2588	ENDDO moist_tile_loop_1
2589	!$OMP END PARALLEL DO
2590
2591	!$OMP PARALLEL DO &
2592	!$OMP PRIVATE ( ij )
2593	moist_tile_loop_2: DO ij = 1 , grid%num_tiles
2594
2595	CALL wrf_debug ( 200 , ' call rk_update_scalar' )
2596	tenddec = .false.
2597
2598	BENCH_START(update_scal_tim)
2599	CALL rk_update_scalar( scs=im, sce=im, &
2600	scalar_1=moist_old(ims,kms,jms,im), &
2601	scalar_2=moist(ims,kms,jms,im), &
2602	sc_tend=moist_tend(ims,kms,jms,im), &
2603	! advh_t=advh_t(ims,kms,jms,1), &
2604	! advz_t=advz_t(ims,kms,jms,1), &
2605	advect_tend=advect_tend, &
2606	h_tendency=h_tendency, z_tendency=z_tendency, &
2607	msftx=grid%msftx,msfty=grid%msfty, &
2608	mu_old=grid%mu_1, mu_new=grid%mu_2, mu_base=grid%mub, &
2609	rk_step=rk_step, dt=dt_rk, spec_zone=grid%spec_zone, &
2610	config_flags=config_flags, tenddec=tenddec, &
2611	ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, &
2612	ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, &
2613	its=grid%i_start(ij), ite=grid%i_end(ij), &
2614	jts=grid%j_start(ij), jte=grid%j_end(ij), &
2615	kts=k_start , kte=k_end )
2616	BENCH_END(update_scal_tim)
2617	PRINT *,' Lluis after rk_update_scalar moist im=',im, &
2618	' moist_old moist_tend moist________'
2619	DO iz = kms, kme
2620	PRINT *,moist_old(config_flags%i_check_point, iz, &
2621	config_flags%j_check_point,im), &
2622	moist_tend(config_flags%i_check_point, iz, &
2623	config_flags%j_check_point,im), &
2624	moist(config_flags%i_check_point, iz, &
2625	config_flags%j_check_point,im)
2626	END DO
2627
2628	BENCH_START(flow_depbdy_tim)
2629	IF( config_flags%specified ) THEN
2630	IF(im .ne. P_QV)THEN
2631	CALL flow_dep_bdy ( moist(ims,kms,jms,im), &
2632	grid%ru_m, grid%rv_m, config_flags, &
2633	grid%spec_zone, &
2634	ids,ide, jds,jde, kds,kde, &
2635	ims,ime, jms,jme, kms,kme, &
2636	ips,ipe, jps,jpe, kps,kpe, &
2637	grid%i_start(ij), grid%i_end(ij), &
2638	grid%j_start(ij), grid%j_end(ij), &
2639	k_start, k_end )
2640	ENDIF
2641	ENDIF
2642	BENCH_END(flow_depbdy_tim)
2643
2644	ENDDO moist_tile_loop_2
2645	!$OMP END PARALLEL DO
2646
2647	ENDIF !-- if (grid%adv_moist_cond .or. im==p_qv ) then
2648
2649	ENDDO moist_variable_loop
2650
2651	ENDIF moist_scalar_advance
2652
2653	BENCH_START(tke_adv_tim)
2654	TKE_advance: IF (config_flags%km_opt .eq. 2) then
2655	#ifdef DM_PARALLEL
2656	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
2657	# include "HALO_EM_TKE_ADVECT_3.inc"
2658	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
2659	# include "HALO_EM_TKE_ADVECT_5.inc"
2660	ELSE
2661	WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order
2662	CALL wrf_error_fatal(TRIM(wrf_err_message))
2663	ENDIF
2664	#endif
2665	!$OMP PARALLEL DO &
2666	!$OMP PRIVATE ( ij )
2667	tke_tile_loop_1: DO ij = 1 , grid%num_tiles
2668
2669	CALL wrf_debug ( 200 , ' call rk_scalar_tend for tke' )
2670	tenddec = .false.
2671	CALL rk_scalar_tend ( 1, 1, config_flags, tenddec, &
2672	rk_step, dt_rk, &
2673	grid%ru_m, grid%rv_m, grid%ww_m, &
2674	grid%muts, grid%mub, grid%mu_1, &
2675	grid%alt, &
2676	grid%tke_1, &
2677	grid%tke_2, &
2678	tke_tend(ims,kms,jms), &
2679	advect_tend,h_tendency,z_tendency,grid%rqvften, &
2680	grid%qv_base, .false., grid%fnm, grid%fnp, &
2681	grid%msfux,grid%msfuy, grid%msfvx, grid%msfvx_inv, &
2682	grid%msfvy, grid%msftx,grid%msfty, &
2683	grid%rdx, grid%rdy, grid%rdn, grid%rdnw, grid%khdif, &
2684	grid%kvdif, grid%xkhh, &
2685	grid%diff_6th_opt, grid%diff_6th_factor, &
2686	config_flags%tke_adv_opt, &
2687	ids, ide, jds, jde, kds, kde, &
2688	ims, ime, jms, jme, kms, kme, &
2689	grid%i_start(ij), grid%i_end(ij), &
2690	grid%j_start(ij), grid%j_end(ij), &
2691	k_start , k_end )
2692
2693	ENDDO tke_tile_loop_1
2694	!$OMP END PARALLEL DO
2695
2696	!$OMP PARALLEL DO &
2697	!$OMP PRIVATE ( ij )
2698	tke_tile_loop_2: DO ij = 1 , grid%num_tiles
2699
2700	CALL wrf_debug ( 200 , ' call rk_update_scalar' )
2701	tenddec = .false.
2702	CALL rk_update_scalar( scs=1, sce=1, &
2703	scalar_1=grid%tke_1, &
2704	scalar_2=grid%tke_2, &
2705	sc_tend=tke_tend(ims,kms,jms), &
2706	! advh_t=advh_t(ims,kms,jms,1), &
2707	! advz_t=advz_t(ims,kms,jms,1), &
2708	advect_tend=advect_tend, &
2709	h_tendency=h_tendency, z_tendency=z_tendency, &
2710	msftx=grid%msftx,msfty=grid%msfty, &
2711	mu_old=grid%mu_1, mu_new=grid%mu_2, mu_base=grid%mub, &
2712	rk_step=rk_step, dt=dt_rk, spec_zone=grid%spec_zone, &
2713	config_flags=config_flags, tenddec=tenddec, &
2714	ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, &
2715	ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, &
2716	its=grid%i_start(ij), ite=grid%i_end(ij), &
2717	jts=grid%j_start(ij), jte=grid%j_end(ij), &
2718	kts=k_start , kte=k_end )
2719
2720	! bound the tke (greater than 0, less than tke_upper_bound)
2721
2722	CALL bound_tke( grid%tke_2, grid%tke_upper_bound, &
2723	ids, ide, jds, jde, kds, kde, &
2724	ims, ime, jms, jme, kms, kme, &
2725	grid%i_start(ij), grid%i_end(ij), &
2726	grid%j_start(ij), grid%j_end(ij), &
2727	k_start , k_end )
2728
2729	IF( config_flags%specified .or. config_flags%nested ) THEN
2730	CALL flow_dep_bdy ( grid%tke_2, &
2731	grid%ru_m, grid%rv_m, config_flags, &
2732	grid%spec_zone, &
2733	ids,ide, jds,jde, kds,kde, & ! domain dims
2734	ims,ime, jms,jme, kms,kme, & ! memory dims
2735	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
2736	grid%i_start(ij), grid%i_end(ij), &
2737	grid%j_start(ij), grid%j_end(ij), &
2738	k_start, k_end )
2739	ENDIF
2740	ENDDO tke_tile_loop_2
2741	!$OMP END PARALLEL DO
2742
2743	ENDIF TKE_advance
2744	BENCH_END(tke_adv_tim)
2745
2746	#ifdef WRF_CHEM
2747	! next the chemical species
2748	BENCH_START(chem_adv_tim)
2749	chem_scalar_advance: IF (num_3d_c >= PARAM_FIRST_SCALAR) THEN
2750
2751	chem_variable_loop: DO ic = PARAM_FIRST_SCALAR, num_3d_c
2752
2753	!$OMP PARALLEL DO &
2754	!$OMP PRIVATE ( ij )
2755	chem_tile_loop_1: DO ij = 1 , grid%num_tiles
2756
2757	CALL wrf_debug ( 200 , ' call rk_scalar_tend in chem_tile_loop_1' )
2758	tenddec = (( config_flags%chemdiag == USECHEMDIAG ) .and. &
2759	( adv_ct_indices(ic) >= PARAM_FIRST_SCALAR ))
2760	CALL rk_scalar_tend ( ic, ic, config_flags, tenddec, &
2761	rk_step, dt_rk, &
2762	grid%ru_m, grid%rv_m, grid%ww_m, &
2763	grid%muts, grid%mub, grid%mu_1, &
2764	grid%alt, &
2765	chem_old(ims,kms,jms,ic), &
2766	chem(ims,kms,jms,ic), &
2767	chem_tend(ims,kms,jms,ic), &
2768	advect_tend,h_tendency,z_tendency,grid%rqvften, &
2769	grid%qv_base, .false., grid%fnm, grid%fnp, &
2770	grid%msfux,grid%msfuy, grid%msfvx, grid%msfvx_inv, &
2771	grid%msfvy, grid%msftx,grid%msfty, &
2772	grid%rdx, grid%rdy, grid%rdn, grid%rdnw, &
2773	grid%khdif, grid%kvdif, grid%xkhh, &
2774	grid%diff_6th_opt, grid%diff_6th_factor, &
2775	config_flags%chem_adv_opt, &
2776	ids, ide, jds, jde, kds, kde, &
2777	ims, ime, jms, jme, kms, kme, &
2778	grid%i_start(ij), grid%i_end(ij), &
2779	grid%j_start(ij), grid%j_end(ij), &
2780	k_start , k_end )
2781	!
2782	! Currently, chemistry species with specified boundaries (i.e. the mother
2783	! domain) are being over written by flow_dep_bdy_chem. So, relax_bdy and
2784	! spec_bdy are only called for nests. For boundary conditions from global model or larger domain,
2785	! chem is uncoupled, and only used for one row/column on inflow (if have_bcs_chem=.true.)
2786	!
2787	IF( ( config_flags%nested ) .and. rk_step == 1 ) THEN
2788	IF(ic.eq.1)CALL wrf_debug ( 10 , ' have_bcs_chem' )
2789	CALL relax_bdy_scalar ( chem_tend(ims,kms,jms,ic), &
2790	chem(ims,kms,jms,ic), grid%mut, &
2791	chem_bxs(jms,kms,1,ic),chem_bxe(jms,kms,1,ic), &
2792	chem_bys(ims,kms,1,ic),chem_bye(ims,kms,1,ic), &
2793	chem_btxs(jms,kms,1,ic),chem_btxe(jms,kms,1,ic), &
2794	chem_btys(ims,kms,1,ic),chem_btye(ims,kms,1,ic), &
2795	config_flags%spec_bdy_width, grid%spec_zone, grid%relax_zone, &
2796	grid%dtbc, grid%fcx, grid%gcx, &
2797	config_flags, &
2798	ids,ide, jds,jde, kds,kde, &
2799	ims,ime, jms,jme, kms,kme, &
2800	ips,ipe, jps,jpe, kps,kpe, &
2801	grid%i_start(ij), grid%i_end(ij), &
2802	grid%j_start(ij), grid%j_end(ij), &
2803	k_start, k_end )
2804	CALL spec_bdy_scalar ( chem_tend(ims,kms,jms,ic), &
2805	chem_bxs(jms,kms,1,ic),chem_bxe(jms,kms,1,ic), &
2806	chem_bys(ims,kms,1,ic),chem_bye(ims,kms,1,ic), &
2807	chem_btxs(jms,kms,1,ic),chem_btxe(jms,kms,1,ic), &
2808	chem_btys(ims,kms,1,ic),chem_btye(ims,kms,1,ic), &
2809	config_flags%spec_bdy_width, grid%spec_zone, &
2810	config_flags, &
2811	ids,ide, jds,jde, kds,kde, &
2812	ims,ime, jms,jme, kms,kme, &
2813	ips,ipe, jps,jpe, kps,kpe, &
2814	grid%i_start(ij), grid%i_end(ij), &
2815	grid%j_start(ij), grid%j_end(ij), &
2816	k_start, k_end )
2817	ENDIF
2818
2819	ENDDO chem_tile_loop_1
2820	!$OMP END PARALLEL DO
2821
2822	!$OMP PARALLEL DO &
2823	!$OMP PRIVATE ( ij )
2824
2825	chem_tile_loop_2: DO ij = 1 , grid%num_tiles
2826
2827	CALL wrf_debug ( 200 , ' call rk_update_scalar' )
2828	tenddec = (( config_flags%chemdiag == USECHEMDIAG ) .and. &
2829	( adv_ct_indices(ic) >= PARAM_FIRST_SCALAR ))
2830	CALL rk_update_scalar( scs=ic, sce=ic, &
2831	scalar_1=chem_old(ims,kms,jms,ic), &
2832	scalar_2=chem(ims,kms,jms,ic), &
2833	sc_tend=chem_tend(ims,kms,jms,ic), &
2834	advh_t=advh_ct(ims,kms,jms,adv_ct_indices(ic)), &
2835	advz_t=advz_ct(ims,kms,jms,adv_ct_indices(ic)), &
2836	advect_tend=advect_tend, &
2837	h_tendency=h_tendency, z_tendency=z_tendency, &
2838	msftx=grid%msftx,msfty=grid%msfty, &
2839	mu_old=grid%mu_1, mu_new=grid%mu_2, mu_base=grid%mub, &
2840	rk_step=rk_step, dt=dt_rk, spec_zone=grid%spec_zone, &
2841	config_flags=config_flags, tenddec=tenddec, &
2842	ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, &
2843	ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, &
2844	its=grid%i_start(ij), ite=grid%i_end(ij), &
2845	jts=grid%j_start(ij), jte=grid%j_end(ij), &
2846	kts=k_start , kte=k_end )
2847
2848	IF( config_flags%specified ) THEN
2849	CALL flow_dep_bdy_chem( chem(ims,kms,jms,ic), &
2850	chem_bxs(jms,kms,1,ic), chem_btxs(jms,kms,1,ic), &
2851	chem_bxe(jms,kms,1,ic), chem_btxe(jms,kms,1,ic), &
2852	chem_bys(ims,kms,1,ic), chem_btys(ims,kms,1,ic), &
2853	chem_bye(ims,kms,1,ic), chem_btye(ims,kms,1,ic), &
2854	dt_rk+grid%dtbc, &
2855	config_flags%spec_bdy_width,grid%z, &
2856	grid%have_bcs_chem, &
2857	grid%ru_m, grid%rv_m, config_flags,grid%alt, &
2858	grid%t_1,grid%pb,grid%p,t0,p1000mb,rcp,grid%ph_2,grid%phb,g, &
2859	grid%spec_zone,ic, &
2860	ids,ide, jds,jde, kds,kde, & ! domain dims
2861	ims,ime, jms,jme, kms,kme, & ! memory dims
2862	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
2863	grid%i_start(ij), grid%i_end(ij), &
2864	grid%j_start(ij), grid%j_end(ij), &
2865	k_start, k_end )
2866	ENDIF
2867	ENDDO chem_tile_loop_2
2868	!$OMP END PARALLEL DO
2869
2870	ENDDO chem_variable_loop
2871	ENDIF chem_scalar_advance
2872	BENCH_END(chem_adv_tim)
2873	#endif
2874	! next the chemical species
2875	BENCH_START(tracer_adv_tim)
2876	tracer_advance: IF (num_tracer >= PARAM_FIRST_SCALAR) THEN
2877
2878	tracer_variable_loop: DO ic = PARAM_FIRST_SCALAR, num_tracer
2879
2880	!$OMP PARALLEL DO &
2881	!$OMP PRIVATE ( ij )
2882	tracer_tile_loop_1: DO ij = 1 , grid%num_tiles
2883
2884	CALL wrf_debug ( 15 , ' call rk_scalar_tend in tracer_tile_loop_1' )
2885	tenddec = .false.
2886	CALL rk_scalar_tend ( ic, ic, config_flags, tenddec, &
2887	rk_step, dt_rk, &
2888	grid%ru_m, grid%rv_m, grid%ww_m, &
2889	grid%muts, grid%mub, grid%mu_1, &
2890	grid%alt, &
2891	tracer_old(ims,kms,jms,ic), &
2892	tracer(ims,kms,jms,ic), &
2893	tracer_tend(ims,kms,jms,ic), &
2894	advect_tend,h_tendency,z_tendency,grid%rqvften, &
2895	grid%qv_base, .false., grid%fnm, grid%fnp, &
2896	grid%msfux,grid%msfuy, grid%msfvx, grid%msfvx_inv, &
2897	grid%msfvy, grid%msftx,grid%msfty, &
2898	grid%rdx, grid%rdy, grid%rdn, grid%rdnw, &
2899	grid%khdif, grid%kvdif, grid%xkhh, &
2900	grid%diff_6th_opt, grid%diff_6th_factor, &
2901	config_flags%tracer_adv_opt, &
2902	ids, ide, jds, jde, kds, kde, &
2903	ims, ime, jms, jme, kms, kme, &
2904	grid%i_start(ij), grid%i_end(ij), &
2905	grid%j_start(ij), grid%j_end(ij), &
2906	k_start , k_end )
2907	!
2908	! Currently, chemistry species with specified boundaries (i.e. the mother
2909	! domain) are being over written by flow_dep_bdy_chem. So, relax_bdy and
2910	! spec_bdy are only called for nests. For boundary conditions from global model or larger domain,
2911	! chem is uncoupled, and only used for one row/column on inflow (if have_bcs_chem=.true.)
2912	!
2913	IF( ( config_flags%nested ) .and. rk_step == 1 ) THEN
2914	IF(ic.eq.1)CALL wrf_debug ( 10 , ' have_bcs_tracer' )
2915	CALL relax_bdy_scalar ( tracer_tend(ims,kms,jms,ic), &
2916	tracer(ims,kms,jms,ic), grid%mut, &
2917	tracer_bxs(jms,kms,1,ic),tracer_bxe(jms,kms,1,ic), &
2918	tracer_bys(ims,kms,1,ic),tracer_bye(ims,kms,1,ic), &
2919	tracer_btxs(jms,kms,1,ic),tracer_btxe(jms,kms,1,ic), &
2920	tracer_btys(ims,kms,1,ic),tracer_btye(ims,kms,1,ic), &
2921	config_flags%spec_bdy_width, grid%spec_zone, grid%relax_zone, &
2922	grid%dtbc, grid%fcx, grid%gcx, &
2923	config_flags, &
2924	ids,ide, jds,jde, kds,kde, &
2925	ims,ime, jms,jme, kms,kme, &
2926	ips,ipe, jps,jpe, kps,kpe, &
2927	grid%i_start(ij), grid%i_end(ij), &
2928	grid%j_start(ij), grid%j_end(ij), &
2929	k_start, k_end )
2930	CALL spec_bdy_scalar ( tracer_tend(ims,kms,jms,ic), &
2931	tracer_bxs(jms,kms,1,ic),tracer_bxe(jms,kms,1,ic), &
2932	tracer_bys(ims,kms,1,ic),tracer_bye(ims,kms,1,ic), &
2933	tracer_btxs(jms,kms,1,ic),tracer_btxe(jms,kms,1,ic), &
2934	tracer_btys(ims,kms,1,ic),tracer_btye(ims,kms,1,ic), &
2935	config_flags%spec_bdy_width, grid%spec_zone, &
2936	config_flags, &
2937	ids,ide, jds,jde, kds,kde, &
2938	ims,ime, jms,jme, kms,kme, &
2939	ips,ipe, jps,jpe, kps,kpe, &
2940	grid%i_start(ij), grid%i_end(ij), &
2941	grid%j_start(ij), grid%j_end(ij), &
2942	k_start, k_end )
2943	ENDIF
2944
2945	ENDDO tracer_tile_loop_1
2946	!$OMP END PARALLEL DO
2947
2948	!$OMP PARALLEL DO &
2949	!$OMP PRIVATE ( ij )
2950
2951	tracer_tile_loop_2: DO ij = 1 , grid%num_tiles
2952
2953	CALL wrf_debug ( 200 , ' call rk_update_scalar' )
2954	tenddec = .false.
2955	CALL rk_update_scalar( scs=ic, sce=ic, &
2956	scalar_1=tracer_old(ims,kms,jms,ic), &
2957	scalar_2=tracer(ims,kms,jms,ic), &
2958	sc_tend=tracer_tend(ims,kms,jms,ic), &
2959	! advh_t=advh_t(ims,kms,jms,1), &
2960	! advz_t=advz_t(ims,kms,jms,1), &
2961	advect_tend=advect_tend, &
2962	h_tendency=h_tendency, z_tendency=z_tendency, &
2963	msftx=grid%msftx,msfty=grid%msfty, &
2964	mu_old=grid%mu_1, mu_new=grid%mu_2, mu_base=grid%mub, &
2965	rk_step=rk_step, dt=dt_rk, spec_zone=grid%spec_zone, &
2966	config_flags=config_flags, tenddec=tenddec, &
2967	ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, &
2968	ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, &
2969	its=grid%i_start(ij), ite=grid%i_end(ij), &
2970	jts=grid%j_start(ij), jte=grid%j_end(ij), &
2971	kts=k_start , kte=k_end )
2972
2973	IF( config_flags%specified ) THEN
2974	#ifdef WRF_CHEM
2975	CALL flow_dep_bdy_tracer( tracer(ims,kms,jms,ic), &
2976	tracer_bxs(jms,kms,1,ic), tracer_btxs(jms,kms,1,ic), &
2977	tracer_bxe(jms,kms,1,ic), tracer_btxe(jms,kms,1,ic), &
2978	tracer_bys(ims,kms,1,ic), tracer_btys(ims,kms,1,ic), &
2979	tracer_bye(ims,kms,1,ic), tracer_btye(ims,kms,1,ic), &
2980	dt_rk+grid%dtbc, &
2981	config_flags%spec_bdy_width,grid%z, &
2982	grid%have_bcs_tracer, &
2983	grid%ru_m, grid%rv_m, config_flags%tracer_opt,grid%alt, &
2984	grid%t_1,grid%pb,grid%p,t0,p1000mb,rcp,grid%ph_2,grid%phb,g, &
2985	grid%spec_zone,ic, &
2986	ids,ide, jds,jde, kds,kde, & ! domain dims
2987	ims,ime, jms,jme, kms,kme, & ! memory dims
2988	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
2989	grid%i_start(ij), grid%i_end(ij), &
2990	grid%j_start(ij), grid%j_end(ij), &
2991	k_start, k_end )
2992	#else
2993	CALL flow_dep_bdy ( tracer(ims,kms,jms,ic), &
2994	grid%ru_m, grid%rv_m, config_flags, &
2995	grid%spec_zone, &
2996	ids,ide, jds,jde, kds,kde, & ! domain dims
2997	ims,ime, jms,jme, kms,kme, & ! memory dims
2998	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
2999	grid%i_start(ij), grid%i_end(ij), &
3000	grid%j_start(ij), grid%j_end(ij), &
3001	k_start, k_end )
3002	#endif
3003	ENDIF
3004	ENDDO tracer_tile_loop_2
3005	!$OMP END PARALLEL DO
3006
3007	ENDDO tracer_variable_loop
3008	ENDIF tracer_advance
3009	BENCH_END(tracer_adv_tim)
3010
3011	! next the other scalar species
3012	other_scalar_advance: IF (num_3d_s >= PARAM_FIRST_SCALAR) THEN
3013
3014	scalar_variable_loop: do is = PARAM_FIRST_SCALAR, num_3d_s
3015	!$OMP PARALLEL DO &
3016	!$OMP PRIVATE ( ij )
3017	scalar_tile_loop_1: DO ij = 1 , grid%num_tiles
3018
3019	CALL wrf_debug ( 200 , ' call rk_scalar_tend' )
3020	tenddec = .false.
3021	CALL rk_scalar_tend ( is, is, config_flags, tenddec, &
3022	rk_step, dt_rk, &
3023	grid%ru_m, grid%rv_m, grid%ww_m, &
3024	grid%muts, grid%mub, grid%mu_1, &
3025	grid%alt, &
3026	scalar_old(ims,kms,jms,is), &
3027	scalar(ims,kms,jms,is), &
3028	scalar_tend(ims,kms,jms,is), &
3029	advect_tend,h_tendency,z_tendency,grid%rqvften, &
3030	grid%qv_base, .false., grid%fnm, grid%fnp, &
3031	grid%msfux,grid%msfuy, grid%msfvx, grid%msfvx_inv, &
3032	grid%msfvy, grid%msftx,grid%msfty, &
3033	grid%rdx, grid%rdy, grid%rdn, grid%rdnw, &
3034	grid%khdif, grid%kvdif, grid%xkhh, &
3035	grid%diff_6th_opt, grid%diff_6th_factor, &
3036	config_flags%scalar_adv_opt, &
3037	ids, ide, jds, jde, kds, kde, &
3038	ims, ime, jms, jme, kms, kme, &
3039	grid%i_start(ij), grid%i_end(ij), &
3040	grid%j_start(ij), grid%j_end(ij), &
3041	k_start , k_end )
3042
3043	IF( config_flags%nested .and. (rk_step == 1) ) THEN
3044
3045	CALL relax_bdy_scalar ( scalar_tend(ims,kms,jms,is), &
3046	scalar(ims,kms,jms,is), grid%mut, &
3047	scalar_bxs(jms,kms,1,is),scalar_bxe(jms,kms,1,is), &
3048	scalar_bys(ims,kms,1,is),scalar_bye(ims,kms,1,is), &
3049	scalar_btxs(jms,kms,1,is),scalar_btxe(jms,kms,1,is), &
3050	scalar_btys(ims,kms,1,is),scalar_btye(ims,kms,1,is), &
3051	config_flags%spec_bdy_width, grid%spec_zone, grid%relax_zone, &
3052	grid%dtbc, grid%fcx, grid%gcx, &
3053	config_flags, &
3054	ids,ide, jds,jde, kds,kde, &
3055	ims,ime, jms,jme, kms,kme, &
3056	ips,ipe, jps,jpe, kps,kpe, &
3057	grid%i_start(ij), grid%i_end(ij), &
3058	grid%j_start(ij), grid%j_end(ij), &
3059	k_start, k_end )
3060
3061	CALL spec_bdy_scalar ( scalar_tend(ims,kms,jms,is), &
3062	scalar_bxs(jms,kms,1,is),scalar_bxe(jms,kms,1,is), &
3063	scalar_bys(ims,kms,1,is),scalar_bye(ims,kms,1,is), &
3064	scalar_btxs(jms,kms,1,is),scalar_btxe(jms,kms,1,is), &
3065	scalar_btys(ims,kms,1,is),scalar_btye(ims,kms,1,is), &
3066	config_flags%spec_bdy_width, grid%spec_zone, &
3067	config_flags, &
3068	ids,ide, jds,jde, kds,kde, &
3069	ims,ime, jms,jme, kms,kme, &
3070	ips,ipe, jps,jpe, kps,kpe, &
3071	grid%i_start(ij), grid%i_end(ij), &
3072	grid%j_start(ij), grid%j_end(ij), &
3073	k_start, k_end )
3074
3075	ENDIF ! b.c test for chem nested boundary condition
3076
3077	ENDDO scalar_tile_loop_1
3078	!$OMP END PARALLEL DO
3079
3080	!$OMP PARALLEL DO &
3081	!$OMP PRIVATE ( ij )
3082	scalar_tile_loop_2: DO ij = 1 , grid%num_tiles
3083
3084	CALL wrf_debug ( 200 , ' call rk_update_scalar' )
3085	tenddec = .false.
3086	CALL rk_update_scalar( scs=is, sce=is, &
3087	scalar_1=scalar_old(ims,kms,jms,is), &
3088	scalar_2=scalar(ims,kms,jms,is), &
3089	sc_tend=scalar_tend(ims,kms,jms,is), &
3090	! advh_t=advh_t(ims,kms,jms,1), &
3091	! advz_t=advz_t(ims,kms,jms,1), &
3092	advect_tend=advect_tend, &
3093	h_tendency=h_tendency, z_tendency=z_tendency, &
3094	msftx=grid%msftx,msfty=grid%msfty, &
3095	mu_old=grid%mu_1, mu_new=grid%mu_2, mu_base=grid%mub, &
3096	rk_step=rk_step, dt=dt_rk, spec_zone=grid%spec_zone, &
3097	config_flags=config_flags, tenddec=tenddec, &
3098	ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, &
3099	ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, &
3100	its=grid%i_start(ij), ite=grid%i_end(ij), &
3101	jts=grid%j_start(ij), jte=grid%j_end(ij), &
3102	kts=k_start , kte=k_end )
3103
3104	IF( config_flags%specified ) THEN
3105
3106	CALL flow_dep_bdy ( scalar(ims,kms,jms,is), &
3107	grid%ru_m, grid%rv_m, config_flags, &
3108	grid%spec_zone, &
3109	ids,ide, jds,jde, kds,kde, & ! domain dims
3110	ims,ime, jms,jme, kms,kme, & ! memory dims
3111	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
3112	grid%i_start(ij), grid%i_end(ij), &
3113	grid%j_start(ij), grid%j_end(ij), &
3114	k_start, k_end )
3115	ENDIF
3116
3117	ENDDO scalar_tile_loop_2
3118	!$OMP END PARALLEL DO
3119
3120	ENDDO scalar_variable_loop
3121
3122	ENDIF other_scalar_advance
3123
3124	! update the pressure and density at the new time level
3125
3126	!$OMP PARALLEL DO &
3127	!$OMP PRIVATE ( ij )
3128	DO ij = 1 , grid%num_tiles
3129
3130	BENCH_START(calc_p_rho_tim)
3131
3132	CALL calc_p_rho_phi( moist, num_3d_m, &
3133	grid%al, grid%alb, grid%mu_2, grid%muts, &
3134	grid%ph_2, grid%p, grid%pb, grid%t_2, &
3135	p0, t0, grid%znu, grid%dnw, grid%rdnw, &
3136	grid%rdn, config_flags%non_hydrostatic, &
3137	ids, ide, jds, jde, kds, kde, &
3138	ims, ime, jms, jme, kms, kme, &
3139	grid%i_start(ij), grid%i_end(ij), &
3140	grid%j_start(ij), grid%j_end(ij), &
3141	k_start , k_end )
3142
3143	BENCH_END(calc_p_rho_tim)
3144
3145	ENDDO
3146	!$OMP END PARALLEL DO
3147
3148	! Reset the boundary conditions if there is another corrector step.
3149	! (rk_step < rk_order), else we'll handle it at the end of everything
3150	! (after the split physics, before exiting the timestep).
3151
3152	rk_step_1_check: IF ( rk_step < rk_order ) THEN
3153
3154	!-----------------------------------------------------------
3155	! rk3 substep polar filter for scalars (moist,chem,scalar)
3156	!-----------------------------------------------------------
3157
3158	IF (config_flags%polar) THEN
3159	IF ( num_3d_m >= PARAM_FIRST_SCALAR ) THEN
3160	CALL wrf_debug ( 200 , ' call filter moist ' )
3161	DO im = PARAM_FIRST_SCALAR, num_3d_m
3162	CALL couple_scalars_for_filter ( FIELD=moist(ims,kms,jms,im) &
3163	,MU=grid%mu_2 , MUB=grid%mub &
3164	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3165	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3166	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe )
3167	CALL pxft ( grid=grid &
3168	,lineno=__LINE__ &
3169	,flag_uv = 0 &
3170	,flag_rurv = 0 &
3171	,flag_wph = 0 &
3172	,flag_ww = 0 &
3173	,flag_t = 0 &
3174	,flag_mu = 0 &
3175	,flag_mut = 0 &
3176	,flag_moist = im &
3177	,flag_chem = 0 &
3178	,flag_scalar = 0 &
3179	,flag_tracer = 0 &
3180	,positive_definite=.FALSE. &
3181	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
3182	,fft_filter_lat = config_flags%fft_filter_lat &
3183	,dclat = dclat &
3184	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3185	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3186	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
3187	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
3188	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
3189	CALL uncouple_scalars_for_filter ( FIELD=moist(ims,kms,jms,im) &
3190	,MU=grid%mu_2 , MUB=grid%mub &
3191	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3192	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3193	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe )
3194	END DO
3195	END IF
3196
3197	IF ( num_3d_c >= PARAM_FIRST_SCALAR ) THEN
3198	CALL wrf_debug ( 200 , ' call filter chem ' )
3199	DO im = PARAM_FIRST_SCALAR, num_3d_c
3200	CALL couple_scalars_for_filter ( FIELD=chem(ims,kms,jms,im) &
3201	,MU=grid%mu_2 , MUB=grid%mub &
3202	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3203	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3204	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe )
3205	CALL pxft ( grid=grid &
3206	,lineno=__LINE__ &
3207	,flag_uv = 0 &
3208	,flag_rurv = 0 &
3209	,flag_wph = 0 &
3210	,flag_ww = 0 &
3211	,flag_t = 0 &
3212	,flag_mu = 0 &
3213	,flag_mut = 0 &
3214	,flag_moist = 0 &
3215	,flag_chem = im &
3216	,flag_tracer = 0 &
3217	,flag_scalar = 0 &
3218	,positive_definite=.FALSE. &
3219	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
3220	,fft_filter_lat = config_flags%fft_filter_lat &
3221	,dclat = dclat &
3222	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3223	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3224	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
3225	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
3226	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
3227	CALL uncouple_scalars_for_filter ( FIELD=chem(ims,kms,jms,im) &
3228	,MU=grid%mu_2 , MUB=grid%mub &
3229	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3230	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3231	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe )
3232	END DO
3233	END IF
3234	IF ( num_tracer >= PARAM_FIRST_SCALAR ) THEN
3235	CALL wrf_debug ( 200 , ' call filter tracer ' )
3236	DO im = PARAM_FIRST_SCALAR, num_tracer
3237	CALL couple_scalars_for_filter ( FIELD=tracer(ims,kms,jms,im) &
3238	,MU=grid%mu_2 , MUB=grid%mub &
3239	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3240	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3241	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe )
3242	CALL pxft ( grid=grid &
3243	,lineno=__LINE__ &
3244	,flag_uv = 0 &
3245	,flag_rurv = 0 &
3246	,flag_wph = 0 &
3247	,flag_ww = 0 &
3248	,flag_t = 0 &
3249	,flag_mu = 0 &
3250	,flag_mut = 0 &
3251	,flag_moist = 0 &
3252	,flag_chem = 0 &
3253	,flag_tracer = im &
3254	,flag_scalar = 0 &
3255	,positive_definite=.FALSE. &
3256	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
3257	,fft_filter_lat = config_flags%fft_filter_lat &
3258	,dclat = dclat &
3259	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3260	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3261	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
3262	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
3263	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
3264	CALL uncouple_scalars_for_filter ( FIELD=tracer(ims,kms,jms,im) &
3265	,MU=grid%mu_2 , MUB=grid%mub &
3266	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3267	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3268	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe )
3269	END DO
3270	END IF
3271
3272	IF ( num_3d_s >= PARAM_FIRST_SCALAR ) THEN
3273	CALL wrf_debug ( 200 , ' call filter scalar ' )
3274	DO im = PARAM_FIRST_SCALAR, num_3d_s
3275	CALL couple_scalars_for_filter ( FIELD=scalar(ims,kms,jms,im) &
3276	,MU=grid%mu_2 , MUB=grid%mub &
3277	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3278	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3279	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe )
3280	CALL pxft ( grid=grid &
3281	,lineno=__LINE__ &
3282	,flag_uv = 0 &
3283	,flag_rurv = 0 &
3284	,flag_wph = 0 &
3285	,flag_ww = 0 &
3286	,flag_t = 0 &
3287	,flag_mu = 0 &
3288	,flag_mut = 0 &
3289	,flag_moist = 0 &
3290	,flag_chem = 0 &
3291	,flag_tracer = 0 &
3292	,flag_scalar = im &
3293	,positive_definite=.FALSE. &
3294	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
3295	,fft_filter_lat = config_flags%fft_filter_lat &
3296	,dclat = dclat &
3297	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3298	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3299	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
3300	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
3301	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
3302	CALL uncouple_scalars_for_filter ( FIELD=scalar(ims,kms,jms,im) &
3303	,MU=grid%mu_2 , MUB=grid%mub &
3304	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3305	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3306	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe )
3307	END DO
3308	END IF
3309	END IF ! polar filter test
3310
3311	!-----------------------------------------------------------
3312	! END rk3 substep polar filter for scalars (moist,chem,scalar)
3313	!-----------------------------------------------------------
3314
3315	!-----------------------------------------------------------
3316	! Stencils for patch communications (WCS, 29 June 2001)
3317	!
3318	! here's where we need a wide comm stencil - these are the
3319	! uncoupled variables so are used for high order calc in
3320	! advection and mixong routines.
3321	!
3322	!
3323	! * * * * * * *
3324	! * * * * * * * * * * * *
3325	! * * * * * * * * * * * * *
3326	! * + * * * + * * * * * + * * *
3327	! * * * * * * * * * * * * *
3328	! * * * * * * * * * * * *
3329	! * * * * * * *
3330	!
3331	! al x
3332	!
3333	! 2D variable
3334	! mu_2 x
3335	!
3336	! (adv order <=4)
3337	! u_2 x
3338	! v_2 x
3339	! w_2 x
3340	! t_2 x
3341	! ph_2 x
3342	!
3343	! (adv order <=6)
3344	! u_2 x
3345	! v_2 x
3346	! w_2 x
3347	! t_2 x
3348	! ph_2 x
3349	!
3350	! 4D variable
3351	! moist x
3352	! chem x
3353	! scalar x
3354
3355	#ifdef DM_PARALLEL
3356	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
3357	# include "HALO_EM_D2_3.inc"
3358	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
3359	# include "HALO_EM_D2_5.inc"
3360	ELSE
3361	WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order
3362	CALL wrf_error_fatal(TRIM(wrf_err_message))
3363	ENDIF
3364	# include "PERIOD_BDY_EM_D.inc"
3365	# include "PERIOD_BDY_EM_MOIST2.inc"
3366	# include "PERIOD_BDY_EM_CHEM2.inc"
3367	# include "PERIOD_BDY_EM_TRACER2.inc"
3368	# include "PERIOD_BDY_EM_SCALAR2.inc"
3369	#endif
3370
3371	BENCH_START(bc_end_tim)
3372	!$OMP PARALLEL DO &
3373	!$OMP PRIVATE ( ij )
3374	tile_bc_loop_1: DO ij = 1 , grid%num_tiles
3375	CALL wrf_debug ( 200 , ' call rk_phys_bc_dry_2' )
3376
3377	CALL rk_phys_bc_dry_2( config_flags, &
3378	grid%u_2, grid%v_2, grid%w_2, &
3379	grid%t_2, grid%ph_2, grid%mu_2, &
3380	ids, ide, jds, jde, kds, kde, &
3381	ims, ime, jms, jme, kms, kme, &
3382	ips, ipe, jps, jpe, kps, kpe, &
3383	grid%i_start(ij), grid%i_end(ij), &
3384	grid%j_start(ij), grid%j_end(ij), &
3385	k_start , k_end )
3386
3387	BENCH_START(diag_w_tim)
3388	IF (.not. config_flags%non_hydrostatic) THEN
3389	CALL diagnose_w( ph_tend, grid%ph_2, grid%ph_1, grid%w_2, grid%muts, dt_rk, &
3390	grid%u_2, grid%v_2, grid%ht, &
3391	grid%cf1, grid%cf2, grid%cf3, grid%rdx, grid%rdy, grid%msftx, grid%msfty, &
3392	ids, ide, jds, jde, kds, kde, &
3393	ims, ime, jms, jme, kms, kme, &
3394	grid%i_start(ij), grid%i_end(ij), &
3395	grid%j_start(ij), grid%j_end(ij), &
3396	k_start , k_end )
3397	ENDIF
3398	BENCH_END(diag_w_tim)
3399
3400	IF (num_3d_m >= PARAM_FIRST_SCALAR) THEN
3401
3402	moisture_loop_bdy_1 : DO im = PARAM_FIRST_SCALAR , num_3d_m
3403
3404	CALL set_physical_bc3d( moist(ims,kms,jms,im), 'p', config_flags, &
3405	ids, ide, jds, jde, kds, kde, &
3406	ims, ime, jms, jme, kms, kme, &
3407	ips, ipe, jps, jpe, kps, kpe, &
3408	grid%i_start(ij), grid%i_end(ij), &
3409	grid%j_start(ij), grid%j_end(ij), &
3410	k_start , k_end )
3411	END DO moisture_loop_bdy_1
3412
3413	ENDIF
3414
3415	IF (num_3d_c >= PARAM_FIRST_SCALAR) THEN
3416
3417	chem_species_bdy_loop_1 : DO ic = PARAM_FIRST_SCALAR , num_3d_c
3418
3419	CALL set_physical_bc3d( chem(ims,kms,jms,ic), 'p', config_flags, &
3420	ids, ide, jds, jde, kds, kde, &
3421	ims, ime, jms, jme, kms, kme, &
3422	ips, ipe, jps, jpe, kps, kpe, &
3423	grid%i_start(ij), grid%i_end(ij), &
3424	grid%j_start(ij), grid%j_end(ij), &
3425	k_start , k_end-1 )
3426
3427	END DO chem_species_bdy_loop_1
3428
3429	END IF
3430
3431	IF (num_tracer >= PARAM_FIRST_SCALAR) THEN
3432
3433	tracer_species_bdy_loop_1 : DO ic = PARAM_FIRST_SCALAR , num_tracer
3434
3435	CALL set_physical_bc3d( tracer(ims,kms,jms,ic), 'p', config_flags, &
3436	ids, ide, jds, jde, kds, kde, &
3437	ims, ime, jms, jme, kms, kme, &
3438	ips, ipe, jps, jpe, kps, kpe, &
3439	grid%i_start(ij), grid%i_end(ij), &
3440	grid%j_start(ij), grid%j_end(ij), &
3441	k_start , k_end-1 )
3442
3443	END DO tracer_species_bdy_loop_1
3444
3445	END IF
3446
3447	IF (num_3d_s >= PARAM_FIRST_SCALAR) THEN
3448
3449	scalar_species_bdy_loop_1 : DO is = PARAM_FIRST_SCALAR , num_3d_s
3450
3451	CALL set_physical_bc3d( scalar(ims,kms,jms,is), 'p', config_flags, &
3452	ids, ide, jds, jde, kds, kde, &
3453	ims, ime, jms, jme, kms, kme, &
3454	ips, ipe, jps, jpe, kps, kpe, &
3455	grid%i_start(ij), grid%i_end(ij), &
3456	grid%j_start(ij), grid%j_end(ij), &
3457	k_start , k_end-1 )
3458
3459	END DO scalar_species_bdy_loop_1
3460
3461	END IF
3462
3463	IF (config_flags%km_opt .eq. 2) THEN
3464
3465	CALL set_physical_bc3d( grid%tke_2 , 'p', config_flags, &
3466	ids, ide, jds, jde, kds, kde, &
3467	ims, ime, jms, jme, kms, kme, &
3468	ips, ipe, jps, jpe, kps, kpe, &
3469	grid%i_start(ij), grid%i_end(ij), &
3470	grid%j_start(ij), grid%j_end(ij), &
3471	k_start , k_end )
3472	END IF
3473
3474	END DO tile_bc_loop_1
3475	!$OMP END PARALLEL DO
3476	BENCH_END(bc_end_tim)
3477
3478
3479	#ifdef DM_PARALLEL
3480
3481	! * * * * *
3482	! * * * * * * * * *
3483	! * + * * + * * * + * *
3484	! * * * * * * * * *
3485	! * * * * *
3486
3487	! moist, chem, scalar, tke x
3488
3489
3490	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
3491	IF ( (config_flags%tke_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3492	# include "HALO_EM_TKE_5.inc"
3493	ELSE
3494	# include "HALO_EM_TKE_3.inc"
3495	ENDIF
3496	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
3497	IF ( (config_flags%tke_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3498	# include "HALO_EM_TKE_7.inc"
3499	ELSE
3500	# include "HALO_EM_TKE_5.inc"
3501	ENDIF
3502	ELSE
3503	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
3504	CALL wrf_error_fatal(TRIM(wrf_err_message))
3505	ENDIF
3506
3507	IF ( num_moist .GE. PARAM_FIRST_SCALAR ) THEN
3508	IF ( config_flags%h_sca_adv_order <= 4 ) THEN
3509	IF ( (config_flags%moist_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3510	# include "HALO_EM_MOIST_E_5.inc"
3511	ELSE
3512	# include "HALO_EM_MOIST_E_3.inc"
3513	END IF
3514	ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN
3515	IF ( (config_flags%moist_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3516	# include "HALO_EM_MOIST_E_7.inc"
3517	ELSE
3518	# include "HALO_EM_MOIST_E_5.inc"
3519	END IF
3520	ELSE
3521	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
3522	CALL wrf_error_fatal(TRIM(wrf_err_message))
3523	ENDIF
3524	ENDIF
3525	IF ( num_chem >= PARAM_FIRST_SCALAR ) THEN
3526	IF ( config_flags%h_sca_adv_order <= 4 ) THEN
3527	IF ( (config_flags%chem_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3528	# include "HALO_EM_CHEM_E_5.inc"
3529	ELSE
3530	# include "HALO_EM_CHEM_E_3.inc"
3531	ENDIF
3532	ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN
3533	IF ( (config_flags%chem_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3534	# include "HALO_EM_CHEM_E_7.inc"
3535	ELSE
3536	# include "HALO_EM_CHEM_E_5.inc"
3537	ENDIF
3538	ELSE
3539	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
3540	CALL wrf_error_fatal(TRIM(wrf_err_message))
3541	ENDIF
3542	ENDIF
3543	IF ( num_tracer >= PARAM_FIRST_SCALAR ) THEN
3544	IF ( config_flags%h_sca_adv_order <= 4 ) THEN
3545	IF ( (config_flags%tracer_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3546	# include "HALO_EM_TRACER_E_5.inc"
3547	ELSE
3548	# include "HALO_EM_TRACER_E_3.inc"
3549	ENDIF
3550	ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN
3551	IF ( (config_flags%tracer_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3552	# include "HALO_EM_TRACER_E_7.inc"
3553	ELSE
3554	# include "HALO_EM_TRACER_E_5.inc"
3555	ENDIF
3556	ELSE
3557	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
3558	CALL wrf_error_fatal(TRIM(wrf_err_message))
3559	ENDIF
3560	ENDIF
3561	IF ( num_scalar >= PARAM_FIRST_SCALAR ) THEN
3562	IF ( config_flags%h_sca_adv_order <= 4 ) THEN
3563	IF ( (config_flags%scalar_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3564	# include "HALO_EM_SCALAR_E_5.inc"
3565	ELSE
3566	# include "HALO_EM_SCALAR_E_3.inc"
3567	ENDIF
3568	ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN
3569	IF ( (config_flags%scalar_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3570	# include "HALO_EM_SCALAR_E_7.inc"
3571	ELSE
3572	# include "HALO_EM_SCALAR_E_5.inc"
3573	ENDIF
3574	ELSE
3575	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
3576	CALL wrf_error_fatal(TRIM(wrf_err_message))
3577	ENDIF
3578	ENDIF
3579	#endif
3580
3581	ENDIF rk_step_1_check
3582
3583
3584	!**********************************************************
3585	!
3586	! end of RK predictor-corrector loop
3587	!
3588	!**********************************************************
3589
3590	END DO Runge_Kutta_loop
3591	#ifdef LMDZ1
3592	WRITE(message, *)' dyn_em: after runge_kutta_loop'
3593	CALL wrf_debug(200, message)
3594	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
3595	' u_tend: ', ru_tendf(im2,1,jm2)
3596	CALL wrf_debug(200, message)
3597	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
3598	'p sfc: ',p8w(im2,kms,jm2)
3599	CALL wrf_debug(200, message)
3600	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
3601	CALL wrf_debug(200, message)
3602	#endif
3603
3604	IF (config_flags%do_avgflx_em .EQ. 1) THEN
3605	! Reinitialize time-averaged fluxes if history output was written after the previous time step:
3606	CALL WRFU_ALARMGET(grid%alarms( HISTORY_ALARM ),prevringtime=temp_time)
3607	CALL domain_clock_get ( grid, current_time=CurrTime, &
3608	current_timestr=message2 )
3609	! use overloaded -, .LT. operator to check whether to initialize avgflx:
3610	! reinitialize after each history output (detect this here by comparing current time
3611	! against last history time and time step - this code follows what's done in adapt_timestep_em):
3612	WRITE ( message , FMT = '("solve_em: old_dt =",g15.6,", dt=",g15.6," on domain ",I3)' ) &
3613	& old_dt,grid%dt,grid%id
3614	CALL wrf_debug(200,message)
3615	old_dt=min(old_dt,grid%dt)
3616	num = INT(old_dt * precision)
3617	den = precision
3618	CALL WRFU_TimeIntervalSet(dtInterval, Sn=num, Sd=den)
3619	IF (CurrTime .lt. temp_time + dtInterval) THEN
3620	WRITE ( message , FMT = '("solve_em: initializing avgflx at time ",A," on domain ",I3)' ) &
3621	& TRIM(message2), grid%id
3622	CALL wrf_message(trim(message))
3623	grid%avgflx_count = 0
3624	!tile-loop for zero_avgflx
3625	!$OMP PARALLEL DO &
3626	!$OMP PRIVATE ( ij )
3627
3628	DO ij = 1 , grid%num_tiles
3629	CALL wrf_debug(200,'In solve_em, before zero_avgflx call')
3630	CALL zero_avgflx(grid%avgflx_rum,grid%avgflx_rvm,grid%avgflx_wwm, &
3631	& ids, ide, jds, jde, kds, kde, &
3632	& ims, ime, jms, jme, kms, kme, &
3633	& grid%i_start(ij), grid%i_end(ij), grid%j_start(ij), grid%j_end(ij), &
3634	& k_start , k_end, f_flux, &
3635	& grid%avgflx_cfu1,grid%avgflx_cfd1,grid%avgflx_dfu1, &
3636	& grid%avgflx_efu1,grid%avgflx_dfd1,grid%avgflx_efd1 )
3637	CALL wrf_debug(200,'In solve_em, after zero_avgflx call')
3638	ENDDO
3639	ENDIF
3640
3641	! Update avgflx quantities
3642	!tile-loop for upd_avgflx
3643	!$OMP PARALLEL DO &
3644	!$OMP PRIVATE ( ij )
3645
3646	DO ij = 1 , grid%num_tiles
3647	CALL wrf_debug(200,'In solve_em, before upd_avgflx call')
3648	CALL upd_avgflx(grid%avgflx_count,grid%avgflx_rum,grid%avgflx_rvm,grid%avgflx_wwm, &
3649	& grid%ru_m, grid%rv_m, grid%ww_m, &
3650	& ids, ide, jds, jde, kds, kde, &
3651	& ims, ime, jms, jme, kms, kme, &
3652	& grid%i_start(ij), grid%i_end(ij), grid%j_start(ij), grid%j_end(ij), &
3653	& k_start , k_end, f_flux, &
3654	& grid%cfu1,grid%cfd1,grid%dfu1,grid%efu1,grid%dfd1,grid%efd1, &
3655	& grid%avgflx_cfu1,grid%avgflx_cfd1,grid%avgflx_dfu1, &
3656	& grid%avgflx_efu1,grid%avgflx_dfd1,grid%avgflx_efd1 )
3657	CALL wrf_debug(200,'In solve_em, after upd_avgflx call')
3658
3659	ENDDO
3660	grid%avgflx_count = grid%avgflx_count + 1
3661	ENDIF
3662	!
3663	!$OMP PARALLEL DO &
3664	!$OMP PRIVATE ( ij )
3665	DO ij = 1 , grid%num_tiles
3666
3667	BENCH_START(advance_ppt_tim)
3668	CALL wrf_debug ( 200 , ' call advance_ppt' )
3669	CALL advance_ppt(grid%rthcuten,grid%rqvcuten,grid%rqccuten,grid%rqrcuten, &
3670	grid%rqicuten,grid%rqscuten, &
3671	grid%rainc,grid%raincv,grid%rainsh,grid%pratec,grid%pratesh, &
3672	grid%nca,grid%htop,grid%hbot,grid%cutop,grid%cubot, &
3673	grid%cuppt, grid%dt, config_flags, &
3674	ids,ide, jds,jde, kds,kde, &
3675	ims,ime, jms,jme, kms,kme, &
3676	grid%i_start(ij), grid%i_end(ij), &
3677	grid%j_start(ij), grid%j_end(ij), &
3678	k_start , k_end )
3679	BENCH_END(advance_ppt_tim)
3680
3681	ENDDO
3682	!$OMP END PARALLEL DO
3683
3684	!<DESCRIPTION>
3685	!<pre>
3686	! (5) time-split physics.
3687	!
3688	! Microphysics are the only time split physics in the WRF model
3689	! at this time. Split-physics begins with the calculation of
3690	! needed diagnostic quantities (pressure, temperature, etc.)
3691	! followed by a call to the microphysics driver,
3692	! and finishes with a clean-up, storing off of a diabatic tendency
3693	! from the moist physics, and a re-calulation of the diagnostic
3694	! quantities pressure and density.
3695	!</pre>
3696	!</DESCRIPTION>
3697
3698	IF( config_flags%specified .or. config_flags%nested ) THEN
3699	sz = grid%spec_zone
3700	ELSE
3701	sz = 0
3702	ENDIF
3703
3704	IF (config_flags%mp_physics /= 0) then
3705
3706	!$OMP PARALLEL DO &
3707	!$OMP PRIVATE ( ij, its, ite, jts, jte )
3708
3709	scalar_tile_loop_1a: DO ij = 1 , grid%num_tiles
3710
3711	IF ( config_flags%periodic_x ) THEN
3712	its = max(grid%i_start(ij),ids)
3713	ite = min(grid%i_end(ij),ide-1)
3714	ELSE
3715	its = max(grid%i_start(ij),ids+sz)
3716	ite = min(grid%i_end(ij),ide-1-sz)
3717	ENDIF
3718	jts = max(grid%j_start(ij),jds+sz)
3719	jte = min(grid%j_end(ij),jde-1-sz)
3720
3721	CALL wrf_debug ( 200 , ' call moist_physics_prep' )
3722	BENCH_START(moist_physics_prep_tim)
3723	CALL moist_physics_prep_em( grid%t_2, grid%t_1, t0, rho, &
3724	grid%al, grid%alb, grid%p, p8w, p0, grid%pb, &
3725	grid%ph_2, grid%phb, th_phy, pi_phy, p_phy, &
3726	grid%z, grid%z_at_w, dz8w, &
3727	dtm, grid%h_diabatic, &
3728	config_flags,grid%fnm, grid%fnp, &
3729	ids, ide, jds, jde, kds, kde, &
3730	ims, ime, jms, jme, kms, kme, &
3731	its, ite, jts, jte, &
3732	k_start , k_end )
3733	BENCH_END(moist_physics_prep_tim)
3734	END DO scalar_tile_loop_1a
3735	!$OMP END PARALLEL DO
3736
3737	CALL wrf_debug ( 200 , ' call microphysics_driver' )
3738
3739	grid%sr = 0.
3740	specified_bdy = config_flags%specified .OR. config_flags%nested
3741	channel_bdy = config_flags%specified .AND. config_flags%periodic_x
3742
3743	BENCH_START(micro_driver_tim)
3744
3745	CALL microphysics_driver( &
3746	& DT=dtm ,DX=grid%dx ,DY=grid%dy &
3747	& ,DZ8W=dz8w ,F_ICE_PHY=grid%f_ice_phy &
3748	& ,ITIMESTEP=grid%itimestep ,LOWLYR=grid%lowlyr &
3749	& ,P8W=p8w ,P=p_phy ,PI_PHY=pi_phy &
3750	& ,RHO=rho ,SPEC_ZONE=grid%spec_zone &
3751	& ,SR=grid%sr ,TH=th_phy &
3752	& ,refl_10cm=grid%refl_10cm & ! hm, 9/22/09 for refl
3753	& ,WARM_RAIN=grid%warm_rain &
3754	& ,T8W=t8w &
3755	& ,CLDFRA=grid%cldfra, EXCH_H=grid%exch_h &
3756	& ,NSOURCE=grid%qndropsource &
3757	#ifdef WRF_CHEM
3758	& ,QLSINK=grid%qlsink,CLDFRA_OLD=grid%cldfra_old &
3759	& ,PRECR=grid%precr, PRECI=grid%preci, PRECS=grid%precs, PRECG=grid%precg &
3760	& ,CHEM_OPT=config_flags%chem_opt, PROGN=config_flags%progn &
3761	#endif
3762	& ,XLAND=grid%xland &
3763	& ,SPECIFIED=specified_bdy, CHANNEL_SWITCH=channel_bdy &
3764	& ,F_RAIN_PHY=grid%f_rain_phy &
3765	& ,F_RIMEF_PHY=grid%f_rimef_phy &
3766	& ,MP_PHYSICS=config_flags%mp_physics &
3767	& ,ID=grid%id &
3768	& ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
3769	& ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
3770	#ifdef RUN_ON_GPU
3771	& ,IPS=ips,IPE=ipe, JPS=jps,JPE=jpe, KPS=kps,KPE=kpe &
3772	#endif
3773	& ,I_START=grid%i_start,I_END=min(grid%i_end, ide-1) &
3774	& ,J_START=grid%j_start,J_END=min(grid%j_end, jde-1) &
3775	& ,KTS=k_start, KTE=min(k_end,kde-1) &
3776	& ,NUM_TILES=grid%num_tiles &
3777	& ,NAER=grid%naer &
3778	! Optional
3779	& , RAINNC=grid%rainnc, RAINNCV=grid%rainncv &
3780	& , SNOWNC=grid%snownc, SNOWNCV=grid%snowncv &
3781	& , GRAUPELNC=grid%graupelnc, GRAUPELNCV=grid%graupelncv & ! for milbrandt2mom
3782	& , HAILNC=grid%hailnc, HAILNCV=grid%hailncv &
3783	& , W=grid%w_2, Z=grid%z, HT=grid%ht &
3784	& , MP_RESTART_STATE=grid%mp_restart_state &
3785	& , TBPVS_STATE=grid%tbpvs_state & ! etampnew
3786	& , TBPVS0_STATE=grid%tbpvs0_state & ! etampnew
3787	& , QV_CURR=moist(ims,kms,jms,P_QV), F_QV=F_QV &
3788	& , QC_CURR=moist(ims,kms,jms,P_QC), F_QC=F_QC &
3789	& , QR_CURR=moist(ims,kms,jms,P_QR), F_QR=F_QR &
3790	& , QI_CURR=moist(ims,kms,jms,P_QI), F_QI=F_QI &
3791	& , QS_CURR=moist(ims,kms,jms,P_QS), F_QS=F_QS &
3792	& , QG_CURR=moist(ims,kms,jms,P_QG), F_QG=F_QG &
3793	& , QH_CURR=moist(ims,kms,jms,P_QH), F_QH=F_QH & ! for milbrandt2mom
3794	& , QNDROP_CURR=scalar(ims,kms,jms,P_QNDROP), F_QNDROP=F_QNDROP &
3795	& , QT_CURR=scalar(ims,kms,jms,P_QT), F_QT=F_QT &
3796	& , QNN_CURR=scalar(ims,kms,jms,P_QNN), F_QNN=F_QNN &
3797	& , QNI_CURR=scalar(ims,kms,jms,P_QNI), F_QNI=F_QNI &
3798	& , QNC_CURR=scalar(ims,kms,jms,P_QNC), F_QNC=F_QNC &
3799	& , QNR_CURR=scalar(ims,kms,jms,P_QNR), F_QNR=F_QNR &
3800	& , QNS_CURR=scalar(ims,kms,jms,P_QNS), F_QNS=F_QNS &
3801	& , QNG_CURR=scalar(ims,kms,jms,P_QNG), F_QNG=F_QNG &
3802	& , QNH_CURR=scalar(ims,kms,jms,P_QNH), F_QNH=F_QNH & ! for milbrandt2mom
3803	! & , QZR_CURR=scalar(ims,kms,jms,P_QZR), F_QZR=F_QZR & ! for milbrandt3mom
3804	! & , QZI_CURR=scalar(ims,kms,jms,P_QZI), F_QZI=F_QZI & ! "
3805	! & , QZS_CURR=scalar(ims,kms,jms,P_QZS), F_QZS=F_QZS & ! "
3806	! & , QZG_CURR=scalar(ims,kms,jms,P_QZG), F_QZG=F_QZG & ! "
3807	! & , QZH_CURR=scalar(ims,kms,jms,P_QZH), F_QZH=F_QZH & ! "
3808	& , qrcuten=grid%rqrcuten, qscuten=grid%rqscuten &
3809	& , qicuten=grid%rqicuten,mu=grid%mut &
3810	& , HAIL=config_flags%gsfcgce_hail & ! for gsfcgce
3811	& , ICE2=config_flags%gsfcgce_2ice & ! for gsfcgce
3812	! & , ccntype=config_flags%milbrandt_ccntype & ! for milbrandt (2mom)
3813	! YLIN
3814	! RI_CURR INPUT
3815	& , RI_CURR=grid%rimi &
3816	)
3817	BENCH_END(micro_driver_tim)
3818
3819	#ifdef LMDZ
3820	grid%h_diabatic = 0.
3821	#endif
3822
3823	#if 0
3824	BENCH_START(microswap_2)
3825	! for load balancing; communication to redistribute the points
3826	IF ( config_flags%mp_physics .EQ. ETAMPNEW ) THEN
3827	#include "SWAP_ETAMP_NEW.inc"
3828	ELSE IF ( config_flags%mp_physics .EQ. WSM3SCHEME ) THEN
3829	#include "SWAP_WSM3.inc"
3830	ENDIF
3831	BENCH_END(microswap_2)
3832	#endif
3833
3834	CALL wrf_debug ( 200 , ' call moist_physics_finish' )
3835	BENCH_START(moist_phys_end_tim)
3836
3837	!$OMP PARALLEL DO &
3838	!$OMP PRIVATE ( ij, its, ite, jts, jte, im, ii, jj, kk )
3839
3840	DO ij = 1 , grid%num_tiles
3841
3842	its = max(grid%i_start(ij),ids)
3843	ite = min(grid%i_end(ij),ide-1)
3844	jts = max(grid%j_start(ij),jds)
3845	jte = min(grid%j_end(ij),jde-1)
3846
3847	CALL microphysics_zero_outb ( &
3848	moist , num_moist , config_flags , &
3849	ids, ide, jds, jde, kds, kde, &
3850	ims, ime, jms, jme, kms, kme, &
3851	its, ite, jts, jte, &
3852	k_start , k_end )
3853
3854	CALL microphysics_zero_outb ( &
3855	scalar , num_scalar , config_flags , &
3856	ids, ide, jds, jde, kds, kde, &
3857	ims, ime, jms, jme, kms, kme, &
3858	its, ite, jts, jte, &
3859	k_start , k_end )
3860
3861	CALL microphysics_zero_outb ( &
3862	chem , num_chem , config_flags , &
3863	ids, ide, jds, jde, kds, kde, &
3864	ims, ime, jms, jme, kms, kme, &
3865	its, ite, jts, jte, &
3866	k_start , k_end )
3867	CALL microphysics_zero_outb ( &
3868	tracer , num_tracer , config_flags , &
3869	ids, ide, jds, jde, kds, kde, &
3870	ims, ime, jms, jme, kms, kme, &
3871	its, ite, jts, jte, &
3872	k_start , k_end )
3873
3874	IF ( config_flags%periodic_x ) THEN
3875	its = max(grid%i_start(ij),ids)
3876	ite = min(grid%i_end(ij),ide-1)
3877	ELSE
3878	its = max(grid%i_start(ij),ids+sz)
3879	ite = min(grid%i_end(ij),ide-1-sz)
3880	ENDIF
3881	jts = max(grid%j_start(ij),jds+sz)
3882	jte = min(grid%j_end(ij),jde-1-sz)
3883
3884	CALL microphysics_zero_outa ( &
3885	moist , num_moist , config_flags , &
3886	ids, ide, jds, jde, kds, kde, &
3887	ims, ime, jms, jme, kms, kme, &
3888	its, ite, jts, jte, &
3889	k_start , k_end )
3890
3891	CALL microphysics_zero_outa ( &
3892	scalar , num_scalar , config_flags , &
3893	ids, ide, jds, jde, kds, kde, &
3894	ims, ime, jms, jme, kms, kme, &
3895	its, ite, jts, jte, &
3896	k_start , k_end )
3897
3898	CALL microphysics_zero_outa ( &
3899	chem , num_chem , config_flags , &
3900	ids, ide, jds, jde, kds, kde, &
3901	ims, ime, jms, jme, kms, kme, &
3902	its, ite, jts, jte, &
3903	k_start , k_end )
3904
3905	CALL microphysics_zero_outa ( &
3906	tracer , num_tracer , config_flags , &
3907	ids, ide, jds, jde, kds, kde, &
3908	ims, ime, jms, jme, kms, kme, &
3909	its, ite, jts, jte, &
3910	k_start , k_end )
3911
3912	CALL moist_physics_finish_em( grid%t_2, grid%t_1, t0, grid%muts, th_phy, &
3913	grid%h_diabatic, dtm, config_flags, &
3914	#if ( WRF_DFI_RADAR == 1 )
3915	grid%dfi_tten_rad,grid%dfi_stage, &
3916	#endif
3917	ids, ide, jds, jde, kds, kde, &
3918	ims, ime, jms, jme, kms, kme, &
3919	its, ite, jts, jte, &
3920	k_start , k_end )
3921
3922	END DO
3923	!$OMP END PARALLEL DO
3924
3925	ENDIF ! microphysics test
3926
3927	!-----------------------------------------------------------
3928	! filter for moist variables post-microphysics and end of timestep
3929	!-----------------------------------------------------------
3930
3931	IF (config_flags%polar) THEN
3932	IF ( num_3d_m >= PARAM_FIRST_SCALAR ) THEN
3933	CALL wrf_debug ( 200 , ' call filter moist' )
3934	DO im = PARAM_FIRST_SCALAR, num_3d_m
3935	DO jj = jps, MIN(jpe,jde-1)
3936	DO kk = kps, MIN(kpe,kde-1)
3937	DO ii = ips, MIN(ipe,ide-1)
3938	moist(ii,kk,jj,im)=moist(ii,kk,jj,im)*(grid%mu_2(ii,jj)+grid%mub(ii,jj))
3939	ENDDO
3940	ENDDO
3941	ENDDO
3942
3943	CALL pxft ( grid=grid &
3944	,lineno=__LINE__ &
3945	,flag_uv = 0 &
3946	,flag_rurv = 0 &
3947	,flag_wph = 0 &
3948	,flag_ww = 0 &
3949	,flag_t = 0 &
3950	,flag_mu = 0 &
3951	,flag_mut = 0 &
3952	,flag_moist = im &
3953	,flag_chem = 0 &
3954	,flag_tracer = 0 &
3955	,flag_scalar = 0 &
3956	,positive_definite=.FALSE. &
3957	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
3958	,fft_filter_lat = config_flags%fft_filter_lat &
3959	,dclat = dclat &
3960	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3961	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3962	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
3963	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
3964	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
3965
3966	DO jj = jps, MIN(jpe,jde-1)
3967	DO kk = kps, MIN(kpe,kde-1)
3968	DO ii = ips, MIN(ipe,ide-1)
3969	moist(ii,kk,jj,im)=moist(ii,kk,jj,im)/(grid%mu_2(ii,jj)+grid%mub(ii,jj))
3970	ENDDO
3971	ENDDO
3972	ENDDO
3973	ENDDO
3974	ENDIF
3975	ENDIF
3976
3977	!-----------------------------------------------------------
3978	! end filter for moist variables post-microphysics and end of timestep
3979	!-----------------------------------------------------------
3980
3981	!$OMP PARALLEL DO &
3982	!$OMP PRIVATE ( ij, its, ite, jts, jte, im, ii, jj, kk )
3983	scalar_tile_loop_1ba: DO ij = 1 , grid%num_tiles
3984
3985	IF ( config_flags%periodic_x ) THEN
3986	its = max(grid%i_start(ij),ids)
3987	ite = min(grid%i_end(ij),ide-1)
3988	ELSE
3989	its = max(grid%i_start(ij),ids+sz)
3990	ite = min(grid%i_end(ij),ide-1-sz)
3991	ENDIF
3992	jts = max(grid%j_start(ij),jds+sz)
3993	jte = min(grid%j_end(ij),jde-1-sz)
3994
3995	CALL calc_p_rho_phi( moist, num_3d_m, &
3996	grid%al, grid%alb, grid%mu_2, grid%muts, &
3997	grid%ph_2, grid%p, grid%pb, grid%t_2, &
3998	p0, t0, grid%znu, grid%dnw, grid%rdnw, &
3999	grid%rdn, config_flags%non_hydrostatic, &
4000	ids, ide, jds, jde, kds, kde, &
4001	ims, ime, jms, jme, kms, kme, &
4002	its, ite, jts, jte, &
4003	k_start , k_end )
4004
4005	END DO scalar_tile_loop_1ba
4006	!$OMP END PARALLEL DO
4007	BENCH_END(moist_phys_end_tim)
4008
4009	IF (.not. config_flags%non_hydrostatic) THEN
4010	#ifdef DM_PARALLEL
4011	# include "HALO_EM_HYDRO_UV.inc"
4012	# include "PERIOD_EM_HYDRO_UV.inc"
4013	#endif
4014	!$OMP PARALLEL DO &
4015	!$OMP PRIVATE ( ij )
4016	DO ij = 1 , grid%num_tiles
4017	CALL diagnose_w( ph_tend, grid%ph_2, grid%ph_1, grid%w_2, grid%muts, dt_rk, &
4018	grid%u_2, grid%v_2, grid%ht, &
4019	grid%cf1, grid%cf2, grid%cf3, grid%rdx, grid%rdy, grid%msftx, grid%msfty, &
4020	ids, ide, jds, jde, kds, kde, &
4021	ims, ime, jms, jme, kms, kme, &
4022	grid%i_start(ij), grid%i_end(ij), &
4023	grid%j_start(ij), grid%j_end(ij), &
4024	k_start , k_end )
4025
4026	END DO
4027	!$OMP END PARALLEL DO
4028
4029	END IF
4030
4031	CALL wrf_debug ( 200 , ' call chem polar filter ' )
4032
4033	!-----------------------------------------------------------
4034	! filter for chem and scalar variables at end of timestep
4035	!-----------------------------------------------------------
4036
4037	IF (config_flags%polar) THEN
4038
4039	IF ( num_3d_c >= PARAM_FIRST_SCALAR ) then
4040	chem_filter_loop: DO im = PARAM_FIRST_SCALAR, num_3d_c
4041	DO jj = jps, MIN(jpe,jde-1)
4042	DO kk = kps, MIN(kpe,kde-1)
4043	DO ii = ips, MIN(ipe,ide-1)
4044	chem(ii,kk,jj,im)=chem(ii,kk,jj,im)*(grid%mu_2(ii,jj)+grid%mub(ii,jj))
4045	ENDDO
4046	ENDDO
4047	ENDDO
4048
4049	CALL pxft ( grid=grid &
4050	,lineno=__LINE__ &
4051	,flag_uv = 0 &
4052	,flag_rurv = 0 &
4053	,flag_wph = 0 &
4054	,flag_ww = 0 &
4055	,flag_t = 0 &
4056	,flag_mu = 0 &
4057	,flag_mut = 0 &
4058	,flag_moist = 0 &
4059	,flag_chem = im &
4060	,flag_tracer = 0 &
4061	,flag_scalar = 0 &
4062	,positive_definite=.FALSE. &
4063	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
4064	,fft_filter_lat = config_flags%fft_filter_lat &
4065	,dclat = dclat &
4066	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
4067	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
4068	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
4069	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
4070	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
4071
4072	DO jj = jps, MIN(jpe,jde-1)
4073	DO kk = kps, MIN(kpe,kde-1)
4074	DO ii = ips, MIN(ipe,ide-1)
4075	chem(ii,kk,jj,im)=chem(ii,kk,jj,im)/(grid%mu_2(ii,jj)+grid%mub(ii,jj))
4076	ENDDO
4077	ENDDO
4078	ENDDO
4079	ENDDO chem_filter_loop
4080	ENDIF
4081	IF ( num_tracer >= PARAM_FIRST_SCALAR ) then
4082	tracer_filter_loop: DO im = PARAM_FIRST_SCALAR, num_tracer
4083	DO jj = jps, MIN(jpe,jde-1)
4084	DO kk = kps, MIN(kpe,kde-1)
4085	DO ii = ips, MIN(ipe,ide-1)
4086	tracer(ii,kk,jj,im)=tracer(ii,kk,jj,im)*(grid%mu_2(ii,jj)+grid%mub(ii,jj))
4087	ENDDO
4088	ENDDO
4089	ENDDO
4090
4091	CALL pxft ( grid=grid &
4092	,lineno=__LINE__ &
4093	,flag_uv = 0 &
4094	,flag_rurv = 0 &
4095	,flag_wph = 0 &
4096	,flag_ww = 0 &
4097	,flag_t = 0 &
4098	,flag_mu = 0 &
4099	,flag_mut = 0 &
4100	,flag_moist = 0 &
4101	,flag_chem = 0 &
4102	,flag_tracer = im &
4103	,flag_scalar = 0 &
4104	,positive_definite=.FALSE. &
4105	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
4106	,fft_filter_lat = config_flags%fft_filter_lat &
4107	,dclat = dclat &
4108	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
4109	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
4110	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
4111	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
4112	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
4113
4114	DO jj = jps, MIN(jpe,jde-1)
4115	DO kk = kps, MIN(kpe,kde-1)
4116	DO ii = ips, MIN(ipe,ide-1)
4117	tracer(ii,kk,jj,im)=tracer(ii,kk,jj,im)/(grid%mu_2(ii,jj)+grid%mub(ii,jj))
4118	ENDDO
4119	ENDDO
4120	ENDDO
4121	ENDDO tracer_filter_loop
4122	ENDIF
4123
4124	IF ( num_3d_s >= PARAM_FIRST_SCALAR ) then
4125	scalar_filter_loop: DO im = PARAM_FIRST_SCALAR, num_3d_s
4126	DO jj = jps, MIN(jpe,jde-1)
4127	DO kk = kps, MIN(kpe,kde-1)
4128	DO ii = ips, MIN(ipe,ide-1)
4129	scalar(ii,kk,jj,im)=scalar(ii,kk,jj,im)*(grid%mu_2(ii,jj)+grid%mub(ii,jj))
4130	ENDDO
4131	ENDDO
4132	ENDDO
4133
4134	CALL pxft ( grid=grid &
4135	,lineno=__LINE__ &
4136	,flag_uv = 0 &
4137	,flag_rurv = 0 &
4138	,flag_wph = 0 &
4139	,flag_ww = 0 &
4140	,flag_t = 0 &
4141	,flag_mu = 0 &
4142	,flag_mut = 0 &
4143	,flag_moist = 0 &
4144	,flag_chem = 0 &
4145	,flag_tracer = 0 &
4146	,flag_scalar = im &
4147	,positive_definite=.FALSE. &
4148	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
4149	,fft_filter_lat = config_flags%fft_filter_lat &
4150	,dclat = dclat &
4151	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
4152	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
4153	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
4154	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
4155	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
4156
4157	DO jj = jps, MIN(jpe,jde-1)
4158	DO kk = kps, MIN(kpe,kde-1)
4159	DO ii = ips, MIN(ipe,ide-1)
4160	scalar(ii,kk,jj,im)=scalar(ii,kk,jj,im)/(grid%mu_2(ii,jj)+grid%mub(ii,jj))
4161	ENDDO
4162	ENDDO
4163	ENDDO
4164	ENDDO scalar_filter_loop
4165	ENDIF
4166	ENDIF
4167
4168	!-----------------------------------------------------------
4169	! end filter for chem and scalar variables at end of timestep
4170	!-----------------------------------------------------------
4171
4172	! We're finished except for boundary condition (and patch) update
4173
4174	! Boundary condition time (or communication time). At this time, we have
4175	! implemented periodic and symmetric physical boundary conditions.
4176
4177	! b.c. routine for data within patch.
4178
4179	! we need to do both time levels of
4180	! data because the time filter only works in the physical solution space.
4181
4182	! First, do patch communications for boundary conditions (periodicity)
4183
4184	!-----------------------------------------------------------
4185	! Stencils for patch communications (WCS, 29 June 2001)
4186	!
4187	! here's where we need a wide comm stencil - these are the
4188	! uncoupled variables so are used for high order calc in
4189	! advection and mixong routines.
4190	!
4191	! * * * * *
4192	! * * * * * * * * *
4193	! * + * * + * * * + * *
4194	! * * * * * * * * *
4195	! * * * * *
4196	!
4197	! grid%u_1 x
4198	! grid%u_2 x
4199	! grid%v_1 x
4200	! grid%v_2 x
4201	! grid%w_1 x
4202	! grid%w_2 x
4203	! grid%t_1 x
4204	! grid%t_2 x
4205	! grid%ph_1 x
4206	! grid%ph_2 x
4207	! grid%tke_1 x
4208	! grid%tke_2 x
4209	!
4210	! 2D variables
4211	! grid%mu_1 x
4212	! grid%mu_2 x
4213	!
4214	! 4D variables
4215	! moist x
4216	! chem x
4217	! scalar x
4218	!----------------------------------------------------------
4219
4220
4221	#ifdef DM_PARALLEL
4222	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
4223	# include "HALO_EM_D3_3.inc"
4224	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
4225	# include "HALO_EM_D3_5.inc"
4226	ELSE
4227	WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order
4228	CALL wrf_error_fatal(TRIM(wrf_err_message))
4229	ENDIF
4230	# include "PERIOD_BDY_EM_D3.inc"
4231	# include "PERIOD_BDY_EM_MOIST.inc"
4232	# include "PERIOD_BDY_EM_CHEM.inc"
4233	# include "PERIOD_BDY_EM_TRACER.inc"
4234	# include "PERIOD_BDY_EM_SCALAR.inc"
4235	#endif
4236
4237	! now set physical b.c on a patch
4238
4239	BENCH_START(bc_2d_tim)
4240	!$OMP PARALLEL DO &
4241	!$OMP PRIVATE ( ij )
4242	tile_bc_loop_2: DO ij = 1 , grid%num_tiles
4243
4244	CALL wrf_debug ( 200 , ' call set_phys_bc_dry_2' )
4245
4246	CALL set_phys_bc_dry_2( config_flags, &
4247	grid%u_1, grid%u_2, grid%v_1, grid%v_2, grid%w_1, grid%w_2, &
4248	grid%t_1, grid%t_2, grid%ph_1, grid%ph_2, grid%mu_1, grid%mu_2, &
4249	ids, ide, jds, jde, kds, kde, &
4250	ims, ime, jms, jme, kms, kme, &
4251	ips, ipe, jps, jpe, kps, kpe, &
4252	grid%i_start(ij), grid%i_end(ij), &
4253	grid%j_start(ij), grid%j_end(ij), &
4254	k_start , k_end )
4255
4256	CALL set_physical_bc3d( grid%tke_1, 'p', config_flags, &
4257	ids, ide, jds, jde, kds, kde, &
4258	ims, ime, jms, jme, kms, kme, &
4259	ips, ipe, jps, jpe, kps, kpe, &
4260	grid%i_start(ij), grid%i_end(ij), &
4261	grid%j_start(ij), grid%j_end(ij), &
4262	k_start , k_end-1 )
4263
4264	CALL set_physical_bc3d( grid%tke_2 , 'p', config_flags, &
4265	ids, ide, jds, jde, kds, kde, &
4266	ims, ime, jms, jme, kms, kme, &
4267	ips, ipe, jps, jpe, kps, kpe, &
4268	grid%i_start(ij), grid%i_end(ij), &
4269	grid%j_start(ij), grid%j_end(ij), &
4270	k_start , k_end )
4271
4272	moisture_loop_bdy_2 : DO im = PARAM_FIRST_SCALAR , num_3d_m
4273
4274	CALL set_physical_bc3d( moist(ims,kms,jms,im), 'p', &
4275	config_flags, &
4276	ids, ide, jds, jde, kds, kde, &
4277	ims, ime, jms, jme, kms, kme, &
4278	ips, ipe, jps, jpe, kps, kpe, &
4279	grid%i_start(ij), grid%i_end(ij), &
4280	grid%j_start(ij), grid%j_end(ij), &
4281	k_start , k_end )
4282
4283	END DO moisture_loop_bdy_2
4284
4285	chem_species_bdy_loop_2 : DO ic = PARAM_FIRST_SCALAR , num_3d_c
4286
4287	CALL set_physical_bc3d( chem(ims,kms,jms,ic) , 'p', config_flags, &
4288	ids, ide, jds, jde, kds, kde, &
4289	ims, ime, jms, jme, kms, kme, &
4290	ips, ipe, jps, jpe, kps, kpe, &
4291	grid%i_start(ij), grid%i_end(ij), &
4292	grid%j_start(ij), grid%j_end(ij), &
4293	k_start , k_end )
4294
4295	END DO chem_species_bdy_loop_2
4296
4297	tracer_species_bdy_loop_2 : DO ic = PARAM_FIRST_SCALAR , num_tracer
4298
4299	CALL set_physical_bc3d( tracer(ims,kms,jms,ic) , 'p', config_flags, &
4300	ids, ide, jds, jde, kds, kde, &
4301	ims, ime, jms, jme, kms, kme, &
4302	ips, ipe, jps, jpe, kps, kpe, &
4303	grid%i_start(ij), grid%i_end(ij), &
4304	grid%j_start(ij), grid%j_end(ij), &
4305	k_start , k_end )
4306
4307	END DO tracer_species_bdy_loop_2
4308
4309	scalar_species_bdy_loop_2 : DO is = PARAM_FIRST_SCALAR , num_3d_s
4310
4311	CALL set_physical_bc3d( scalar(ims,kms,jms,is) , 'p', config_flags, &
4312	ids, ide, jds, jde, kds, kde, &
4313	ims, ime, jms, jme, kms, kme, &
4314	ips, ipe, jps, jpe, kps, kpe, &
4315	grid%i_start(ij), grid%i_end(ij), &
4316	grid%j_start(ij), grid%j_end(ij), &
4317	k_start , k_end )
4318
4319	END DO scalar_species_bdy_loop_2
4320
4321	END DO tile_bc_loop_2
4322	!$OMP END PARALLEL DO
4323	BENCH_END(bc_2d_tim)
4324
4325	IF( config_flags%specified .or. config_flags%nested ) THEN
4326	grid%dtbc = grid%dtbc + grid%dt
4327	ENDIF
4328
4329	! reset surface w for consistency
4330
4331	#ifdef DM_PARALLEL
4332	# include "HALO_EM_C.inc"
4333	# include "PERIOD_BDY_EM_E.inc"
4334	#endif
4335
4336	CALL wrf_debug ( 10 , ' call set_w_surface' )
4337	fill_w_flag = .false.
4338
4339	!$OMP PARALLEL DO &
4340	!$OMP PRIVATE ( ij )
4341	DO ij = 1 , grid%num_tiles
4342	CALL set_w_surface( config_flags, grid%znw, fill_w_flag, &
4343	grid%w_2, grid%ht, grid%u_2, grid%v_2, &
4344	grid%cf1, grid%cf2, grid%cf3, grid%rdx, grid%rdy,&
4345	grid%msftx, grid%msfty, &
4346	ids, ide, jds, jde, kds, kde, &
4347	ims, ime, jms, jme, kms, kme, &
4348	grid%i_start(ij), grid%i_end(ij), &
4349	grid%j_start(ij), grid%j_end(ij), &
4350	k_start, k_end )
4351	! its, ite, jts, jte, k_start, min(k_end,kde-1), &
4352
4353	END DO
4354	!$OMP END PARALLEL DO
4355
4356	! calculate some model diagnostics.
4357
4358	#ifdef LMDZ1
4359	WRITE(message, *)' dyn_em: before diagnostics'
4360	CALL wrf_debug(200, message)
4361	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
4362	' u_tend: ', ru_tendf(im2,1,jm2)
4363	CALL wrf_debug(200, message)
4364	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
4365	'p sfc: ',p8w(im2,kms,jm2)
4366	CALL wrf_debug(200, message)
4367	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
4368	CALL wrf_debug(200, message)
4369	#endif
4370	CALL wrf_debug ( 200 , ' call diagnostic_driver' )
4371
4372	CALL diagnostic_output_calc( &
4373	& DPSDT=grid%dpsdt ,DMUDT=grid%dmudt &
4374	& ,P8W=p8w ,PK1M=grid%pk1m &
4375	& ,MU_2=grid%mu_2 ,MU_2M=grid%mu_2m &
4376	& ,U=grid%u_2 ,V=grid%v_2 &
4377	& ,RAINCV=grid%raincv ,RAINNCV=grid%rainncv &
4378	& ,RAINC=grid%rainc ,RAINNC=grid%rainnc &
4379	& ,I_RAINC=grid%i_rainc ,I_RAINNC=grid%i_rainnc &
4380	& ,HFX=grid%hfx ,SFCEVP=grid%sfcevp ,LH=grid%lh &
4381	& ,DT=grid%dt ,SBW=config_flags%spec_bdy_width &
4382	& ,XTIME=grid%xtime ,T2=grid%t2 &
4383	& ,ACSWUPT=grid%acswupt ,ACSWUPTC=grid%acswuptc &
4384	& ,ACSWDNT=grid%acswdnt ,ACSWDNTC=grid%acswdntc &
4385	& ,ACSWUPB=grid%acswupb ,ACSWUPBC=grid%acswupbc &
4386	& ,ACSWDNB=grid%acswdnb ,ACSWDNBC=grid%acswdnbc &
4387	& ,ACLWUPT=grid%aclwupt ,ACLWUPTC=grid%aclwuptc &
4388	& ,ACLWDNT=grid%aclwdnt ,ACLWDNTC=grid%aclwdntc &
4389	& ,ACLWUPB=grid%aclwupb ,ACLWUPBC=grid%aclwupbc &
4390	& ,ACLWDNB=grid%aclwdnb ,ACLWDNBC=grid%aclwdnbc &
4391	& ,I_ACSWUPT=grid%i_acswupt ,I_ACSWUPTC=grid%i_acswuptc &
4392	& ,I_ACSWDNT=grid%i_acswdnt ,I_ACSWDNTC=grid%i_acswdntc &
4393	& ,I_ACSWUPB=grid%i_acswupb ,I_ACSWUPBC=grid%i_acswupbc &
4394	& ,I_ACSWDNB=grid%i_acswdnb ,I_ACSWDNBC=grid%i_acswdnbc &
4395	& ,I_ACLWUPT=grid%i_aclwupt ,I_ACLWUPTC=grid%i_aclwuptc &
4396	& ,I_ACLWDNT=grid%i_aclwdnt ,I_ACLWDNTC=grid%i_aclwdntc &
4397	& ,I_ACLWUPB=grid%i_aclwupb ,I_ACLWUPBC=grid%i_aclwupbc &
4398	& ,I_ACLWDNB=grid%i_aclwdnb ,I_ACLWDNBC=grid%i_aclwdnbc &
4399	! Selection flag
4400	& ,DIAG_PRINT=config_flags%diag_print &
4401	& ,BUCKET_MM=config_flags%bucket_mm &
4402	& ,BUCKET_J =config_flags%bucket_J &
4403	& ,SNOWNCV=grid%snowncv, SNOW_ACC_NC=grid%snow_acc_nc &
4404	& ,PREC_ACC_C=grid%prec_acc_c &
4405	& ,PREC_ACC_NC=grid%prec_acc_nc &
4406	& ,PREC_ACC_DT=config_flags%prec_acc_dt &
4407	& ,CURR_SECS=curr_secs &
4408	! Dimension arguments
4409	& ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
4410	& ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
4411	& ,IPS=ips,IPE=ipe, JPS=jps,JPE=jpe, KPS=kps,KPE=kpe &
4412	& ,I_START=grid%i_start,I_END=min(grid%i_end, ide-1) &
4413	& ,J_START=grid%j_start,J_END=min(grid%j_end, jde-1) &
4414	& ,KTS=k_start, KTE=min(k_end,kde-1) &
4415	& ,NUM_TILES=grid%num_tiles &
4416	& )
4417
4418	#ifdef DM_PARALLEL
4419	!-----------------------------------------------------------------------
4420	! see above
4421	!--------------------------------------------------------------
4422	CALL wrf_debug ( 200 , ' call HALO_RK_E' )
4423	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
4424	# include "HALO_EM_E_3.inc"
4425	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
4426	# include "HALO_EM_E_5.inc"
4427	ELSE
4428	WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order
4429	CALL wrf_error_fatal(TRIM(wrf_err_message))
4430	ENDIF
4431	#endif
4432
4433	#ifdef DM_PARALLEL
4434	IF ( num_moist >= PARAM_FIRST_SCALAR ) THEN
4435	!-----------------------------------------------------------------------
4436	! see above
4437	!--------------------------------------------------------------
4438	CALL wrf_debug ( 200 , ' call HALO_RK_MOIST' )
4439	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
4440	# include "HALO_EM_MOIST_E_3.inc"
4441	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
4442	# include "HALO_EM_MOIST_E_5.inc"
4443	ELSE
4444	WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order
4445	CALL wrf_error_fatal(TRIM(wrf_err_message))
4446	ENDIF
4447	ENDIF
4448	IF ( num_chem >= PARAM_FIRST_SCALAR ) THEN
4449	!-----------------------------------------------------------------------
4450	! see above
4451	!--------------------------------------------------------------
4452	CALL wrf_debug ( 200 , ' call HALO_RK_CHEM' )
4453	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
4454	# include "HALO_EM_CHEM_E_3.inc"
4455	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
4456	# include "HALO_EM_CHEM_E_5.inc"
4457	ELSE
4458	WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order
4459	CALL wrf_error_fatal(TRIM(wrf_err_message))
4460	ENDIF
4461	ENDIF
4462	IF ( num_tracer >= PARAM_FIRST_SCALAR ) THEN
4463	!-----------------------------------------------------------------------
4464	! see above
4465	!--------------------------------------------------------------
4466	CALL wrf_debug ( 200 , ' call HALO_RK_TRACER' )
4467	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
4468	# include "HALO_EM_TRACER_E_3.inc"
4469	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
4470	# include "HALO_EM_TRACER_E_5.inc"
4471	ELSE
4472	WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order
4473	CALL wrf_error_fatal(TRIM(wrf_err_message))
4474	ENDIF
4475	ENDIF
4476	IF ( num_scalar >= PARAM_FIRST_SCALAR ) THEN
4477	!-----------------------------------------------------------------------
4478	! see above
4479	!--------------------------------------------------------------
4480	CALL wrf_debug ( 200 , ' call HALO_RK_SCALAR' )
4481	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
4482	# include "HALO_EM_SCALAR_E_3.inc"
4483	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
4484	# include "HALO_EM_SCALAR_E_5.inc"
4485	ELSE
4486	WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order
4487	CALL wrf_error_fatal(TRIM(wrf_err_message))
4488	ENDIF
4489	ENDIF
4490	#endif
4491
4492	! Max values of CFL for adaptive time step scheme
4493
4494	DEALLOCATE(max_vert_cfl_tmp)
4495	DEALLOCATE(max_horiz_cfl_tmp)
4496
4497	CALL wrf_debug ( 200 , ' call end of solve_em' )
4498
4499	! Finish timers if compiled with -DBENCH.
4500	#include <bench_solve_em_end.h>
4501
4502	RETURN
4503
4504	END SUBROUTINE solve_em

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

Download in other formats:

Original Format