Context Navigation

solve_em.F @ 354

Last change on this file since 354 was 258, checked in by lfita, 10 years ago
Activating moisture evolutions even for mp_physics=0
File size: 226.7 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_MAPFAC=grid%msft, &
918	!! & WRF_UTEND=grid%ru_tend, WRF_VTEND=grid%rv_tend, &
919	!! & WRF_TTEND=t_tend, &
920	& WRF_UTEND=ru_tendf, WRF_VTEND=rv_tendf, &
921	& WRF_TTEND=t_tendf, &
922	& WRF_MOISTTEND=moist_tend, WRF_PSFCTEND=grid%dpsdt, &
923	& WRF_QVID=P_QV, WRF_QCID=P_QC, WRF_QRID=P_QR, &
924	& WRF_QSID=P_QS, WRF_QIID=P_QI, WRF_QHID=P_QH, WRF_QGID=P_QG, &
925	! L. Fita. July 2013. Now defined as local dummy variables
926	! & WRF_DUDYN=???????????????, WRF_PVTHETA=????????????????????????, &
927	! & WRF_CLESPHY=?????????????, WRF_PRESNIVS=???????????????????????, &
928	& WRF_MAPFT=grid%msft, WRF_MAPFU=grid%msfu, WRF_MAPFV=grid%msfv, &
929	& WRF_DX=grid%dx, WRF_DY=grid%dy, &
930	& WRF_DBG=model_config_rec%debug_level, LANDCAT=mminlu, &
931	& SOILCAT=mminsl, &
932	& WRF_L_PBL=config_flags%lmdz_iflag_pbl, &
933	& WRF_L_CON=config_flags%lmdz_iflag_con, &
934	& WRF_L_THERMALS=config_flags%lmdz_iflag_thermals, &
935	& WRF_L_WAKE=config_flags%lmdz_iflag_wake, &
936	& wrf_nsoillayers=config_flags%num_soil_layers, &
937	& ICHECK_P=config_flags%i_check_point, &
938	& JCHECK_P=config_flags%j_check_point, &
939	& KCHECK_P=config_flags%k_check_point &
940	& )
941	END IF
942	#endif
943
944	#ifdef LMDZ1
945	WRITE(message, *)' dyn_em: post lmdz pre step2 rk_step:', rk_step
946	CALL wrf_debug(200, message)
947	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
948	' u_tend: ', ru_tendf(im2,1,jm2)
949	CALL wrf_debug(200, message)
950	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
951	'p sfc: ',p8w(im2,kms,jm2)
952	CALL wrf_debug(200, message)
953	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
954	CALL wrf_debug(200, message)
955	#endif
956
957	CALL first_rk_step_part2 ( grid, config_flags &
958	, moist , moist_tend &
959	, chem , chem_tend &
960	, tracer, tracer_tend &
961	, scalar , scalar_tend &
962	, fdda3d, fdda2d &
963	, ru_tendf, rv_tendf &
964	, rw_tendf, t_tendf &
965	, ph_tendf, mu_tendf &
966	, tke_tend &
967	, adapt_step_flag , curr_secs &
968	, psim , psih , wspd , gz1oz0 &
969	, br , chklowq &
970	, cu_act_flag , hol , th_phy &
971	, pi_phy , p_phy , grid%t_phy &
972	, u_phy , v_phy &
973	, dz8w , p8w , t8w , rho_phy , rho &
974	, nba_mij, num_nba_mij & !JDM
975	, nba_rij, num_nba_rij & !JDM
976	, ids, ide, jds, jde, kds, kde &
977	, ims, ime, jms, jme, kms, kme &
978	, ips, ipe, jps, jpe, kps, kpe &
979	, imsx, imex, jmsx, jmex, kmsx, kmex &
980	, ipsx, ipex, jpsx, jpex, kpsx, kpex &
981	, imsy, imey, jmsy, jmey, kmsy, kmey &
982	, ipsy, ipey, jpsy, jpey, kpsy, kpey &
983	, k_start , k_end &
984	)
985
986	END IF rk_step_is_one
987
988	#ifdef LMDZ1
989	WRITE(message, *)' dyn_em: post step2 pre rk_tendency rk_step: ',rk_step
990	CALL wrf_debug(200, message)
991	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
992	' u_tend: ', ru_tendf(im2,1,jm2)
993	CALL wrf_debug(200, message)
994	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
995	'p sfc: ',p8w(im2,kms,jm2)
996	CALL wrf_debug(200, message)
997	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
998	CALL wrf_debug(200, message)
999	#endif
1000
1001	BENCH_START(rk_tend_tim)
1002	!$OMP PARALLEL DO &
1003	!$OMP PRIVATE ( ij )
1004	DO ij = 1 , grid%num_tiles
1005
1006	CALL wrf_debug ( 200 , ' call rk_tendency' )
1007	CALL rk_tendency ( config_flags, rk_step &
1008	,grid%ru_tend, grid%rv_tend, rw_tend, ph_tend, t_tend &
1009	,ru_tendf, rv_tendf, rw_tendf, ph_tendf, t_tendf &
1010	,mu_tend, grid%u_save, grid%v_save, w_save, ph_save &
1011	,grid%t_save, mu_save, grid%rthften &
1012	,grid%ru, grid%rv, grid%rw, grid%ww &
1013	,grid%u_2, grid%v_2, grid%w_2, grid%t_2, grid%ph_2 &
1014	,grid%u_1, grid%v_1, grid%w_1, grid%t_1, grid%ph_1 &
1015	,grid%h_diabatic, grid%phb, grid%t_init &
1016	,grid%mu_2, grid%mut, grid%muu, grid%muv, grid%mub &
1017	,grid%al, grid%alt, grid%p, grid%pb, grid%php, cqu, cqv, cqw &
1018	,grid%u_base, grid%v_base, grid%t_base, grid%qv_base, grid%z_base &
1019	,grid%msfux,grid%msfuy, grid%msfvx, grid%msfvx_inv &
1020	,grid%msfvy, grid%msftx,grid%msfty, grid%clat, grid%f, grid%e, grid%sina, grid%cosa &
1021	,grid%fnm, grid%fnp, grid%rdn, grid%rdnw &
1022	,grid%dt, grid%rdx, grid%rdy, grid%khdif, grid%kvdif, grid%xkmh, grid%xkhh &
1023	,grid%diff_6th_opt, grid%diff_6th_factor &
1024	,grid%dampcoef,grid%zdamp,config_flags%damp_opt,config_flags%rad_nudge &
1025	,grid%cf1, grid%cf2, grid%cf3, grid%cfn, grid%cfn1, num_3d_m &
1026	,config_flags%non_hydrostatic, config_flags%top_lid &
1027	,grid%u_frame, grid%v_frame &
1028	,ids, ide, jds, jde, kds, kde &
1029	,ims, ime, jms, jme, kms, kme &
1030	,grid%i_start(ij), grid%i_end(ij) &
1031	,grid%j_start(ij), grid%j_end(ij) &
1032	,k_start, k_end &
1033	,max_vert_cfl_tmp(ij), max_horiz_cfl_tmp(ij) )
1034
1035	END DO
1036	#ifdef LMDZ1
1037	WRITE(message, *)' dyn_em: post rk_tendency'
1038	CALL wrf_debug(200, message)
1039	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1040	' u_tend: ', ru_tendf(im2,1,jm2)
1041	CALL wrf_debug(200, message)
1042	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1043	'p sfc: ',p8w(im2,kms,jm2)
1044	CALL wrf_debug(200, message)
1045	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1046	CALL wrf_debug(200, message)
1047	#endif
1048	!$OMP END PARALLEL DO
1049	BENCH_END(rk_tend_tim)
1050
1051	IF (config_flags%use_adaptive_time_step) THEN
1052	DO ij = 1 , grid%num_tiles
1053	IF (max_horiz_cfl_tmp(ij) .GT. grid%max_horiz_cfl) THEN
1054	grid%max_horiz_cfl = max_horiz_cfl_tmp(ij)
1055	ENDIF
1056	IF (max_vert_cfl_tmp(ij) .GT. grid%max_vert_cfl) THEN
1057	grid%max_vert_cfl = max_vert_cfl_tmp(ij)
1058	ENDIF
1059	END DO
1060
1061	IF (grid%max_horiz_cfl .GT. grid%max_cfl_val) THEN
1062	grid%max_cfl_val = grid%max_horiz_cfl
1063	ENDIF
1064	IF (grid%max_vert_cfl .GT. grid%max_cfl_val) THEN
1065	grid%max_cfl_val = grid%max_vert_cfl
1066	ENDIF
1067	ENDIF
1068
1069	BENCH_START(relax_bdy_dry_tim)
1070	!$OMP PARALLEL DO &
1071	!$OMP PRIVATE ( ij )
1072	DO ij = 1 , grid%num_tiles
1073
1074	IF( (config_flags%specified .or. config_flags%nested) .and. rk_step == 1 ) THEN
1075
1076	CALL relax_bdy_dry ( config_flags, &
1077	grid%u_save, grid%v_save, ph_save, grid%t_save, &
1078	w_save, mu_tend, &
1079	grid%ru, grid%rv, grid%ph_2, grid%t_2, &
1080	grid%w_2, grid%mu_2, grid%mut, &
1081	grid%u_bxs,grid%u_bxe,grid%u_bys,grid%u_bye, &
1082	grid%v_bxs,grid%v_bxe,grid%v_bys,grid%v_bye, &
1083	grid%ph_bxs,grid%ph_bxe,grid%ph_bys,grid%ph_bye, &
1084	grid%t_bxs,grid%t_bxe,grid%t_bys,grid%t_bye, &
1085	grid%w_bxs,grid%w_bxe,grid%w_bys,grid%w_bye, &
1086	grid%mu_bxs,grid%mu_bxe,grid%mu_bys,grid%mu_bye, &
1087	grid%u_btxs,grid%u_btxe,grid%u_btys,grid%u_btye, &
1088	grid%v_btxs,grid%v_btxe,grid%v_btys,grid%v_btye, &
1089	grid%ph_btxs,grid%ph_btxe,grid%ph_btys,grid%ph_btye, &
1090	grid%t_btxs,grid%t_btxe,grid%t_btys,grid%t_btye, &
1091	grid%w_btxs,grid%w_btxe,grid%w_btys,grid%w_btye, &
1092	grid%mu_btxs,grid%mu_btxe,grid%mu_btys,grid%mu_btye, &
1093	config_flags%spec_bdy_width, grid%spec_zone, grid%relax_zone, &
1094	grid%dtbc, grid%fcx, grid%gcx, &
1095	ids,ide, jds,jde, kds,kde, &
1096	ims,ime, jms,jme, kms,kme, &
1097	ips,ipe, jps,jpe, kps,kpe, &
1098	grid%i_start(ij), grid%i_end(ij), &
1099	grid%j_start(ij), grid%j_end(ij), &
1100	k_start, k_end )
1101
1102	ENDIF
1103
1104	CALL rk_addtend_dry( grid%ru_tend, grid%rv_tend, rw_tend, ph_tend, t_tend, &
1105	ru_tendf, rv_tendf, rw_tendf, ph_tendf, t_tendf, &
1106	grid%u_save, grid%v_save, w_save, ph_save, grid%t_save, &
1107	mu_tend, mu_tendf, rk_step, &
1108	grid%h_diabatic, grid%mut, grid%msftx, &
1109	grid%msfty, grid%msfux,grid%msfuy, &
1110	grid%msfvx, grid%msfvx_inv, grid%msfvy, &
1111	ids,ide, jds,jde, kds,kde, &
1112	ims,ime, jms,jme, kms,kme, &
1113	ips,ipe, jps,jpe, kps,kpe, &
1114	grid%i_start(ij), grid%i_end(ij), &
1115	grid%j_start(ij), grid%j_end(ij), &
1116	k_start, k_end )
1117
1118	IF( config_flags%specified .or. config_flags%nested ) THEN
1119	CALL spec_bdy_dry ( config_flags, &
1120	grid%ru_tend, grid%rv_tend, ph_tend, t_tend, &
1121	rw_tend, mu_tend, &
1122	grid%u_bxs,grid%u_bxe,grid%u_bys,grid%u_bye, &
1123	grid%v_bxs,grid%v_bxe,grid%v_bys,grid%v_bye, &
1124	grid%ph_bxs,grid%ph_bxe,grid%ph_bys,grid%ph_bye, &
1125	grid%t_bxs,grid%t_bxe,grid%t_bys,grid%t_bye, &
1126	grid%w_bxs,grid%w_bxe,grid%w_bys,grid%w_bye, &
1127	grid%mu_bxs,grid%mu_bxe,grid%mu_bys,grid%mu_bye, &
1128	grid%u_btxs,grid%u_btxe,grid%u_btys,grid%u_btye, &
1129	grid%v_btxs,grid%v_btxe,grid%v_btys,grid%v_btye, &
1130	grid%ph_btxs,grid%ph_btxe,grid%ph_btys,grid%ph_btye, &
1131	grid%t_btxs,grid%t_btxe,grid%t_btys,grid%t_btye, &
1132	grid%w_btxs,grid%w_btxe,grid%w_btys,grid%w_btye, &
1133	grid%mu_btxs,grid%mu_btxe,grid%mu_btys,grid%mu_btye, &
1134	config_flags%spec_bdy_width, grid%spec_zone, &
1135	ids,ide, jds,jde, kds,kde, & ! domain dims
1136	ims,ime, jms,jme, kms,kme, & ! memory dims
1137	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
1138	grid%i_start(ij), grid%i_end(ij), &
1139	grid%j_start(ij), grid%j_end(ij), &
1140	k_start, k_end )
1141
1142	ENDIF
1143
1144	END DO
1145	!$OMP END PARALLEL DO
1146	BENCH_END(relax_bdy_dry_tim)
1147
1148	!<DESCRIPTION>
1149	!<pre>
1150	! (3) Small (acoustic,sound) steps.
1151	!
1152	! Several acoustic steps are taken each RK pass. A small step
1153	! sequence begins with calculating perturbation variables
1154	! and coupling them to the column dry-air-mass mu
1155	! (call to small_step_prep). This is followed by computing
1156	! coefficients for the vertically implicit part of the
1157	! small timestep (call to calc_coef_w).
1158	!
1159	! The small steps are taken
1160	! in the named loop "small_steps:". In the small_steps loop, first
1161	! the horizontal momentum (u and v) are advanced (call to advance_uv),
1162	! next mu and theta are advanced (call to advance_mu_t) followed by
1163	! advancing w and the geopotential (call to advance_w). Diagnostic
1164	! values for pressure and inverse density are updated at the end of
1165	! each small_step.
1166	!
1167	! The small-step section ends with the change of the perturbation variables
1168	! back to full variables (call to small_step_finish).
1169	!</pre>
1170	!</DESCRIPTION>
1171
1172	BENCH_START(small_step_prep_tim)
1173	!$OMP PARALLEL DO &
1174	!$OMP PRIVATE ( ij )
1175	DO ij = 1 , grid%num_tiles
1176
1177	! Calculate coefficients for the vertically implicit acoustic/gravity wave
1178	! integration. We only need calculate these for the first pass through -
1179	! the predictor step. They are reused as is for the corrector step.
1180	! For third-order RK, we need to recompute these after the first
1181	! predictor because we may have changed the small timestep -> grid%dts.
1182
1183	CALL wrf_debug ( 200 , ' call small_step_prep ' )
1184
1185	#ifdef LMDZ1
1186	WRITE(message, *)' dyn_em: before small_step_prep'
1187	CALL wrf_debug(200, message)
1188	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1189	' u_tend: ', ru_tendf(im2,1,jm2)
1190	CALL wrf_debug(200, message)
1191	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1192	'p sfc: ',p8w(im2,kms,jm2)
1193	CALL wrf_debug(200, message)
1194	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1195	CALL wrf_debug(200, message)
1196	#endif
1197	CALL small_step_prep( grid%u_1,grid%u_2,grid%v_1,grid%v_2,grid%w_1,grid%w_2, &
1198	grid%t_1,grid%t_2,grid%ph_1,grid%ph_2, &
1199	grid%mub, grid%mu_1, grid%mu_2, &
1200	grid%muu, muus, grid%muv, muvs, &
1201	grid%mut, grid%muts, grid%mudf, &
1202	grid%u_save, grid%v_save, w_save, &
1203	grid%t_save, ph_save, mu_save, &
1204	grid%ww, ww1, &
1205	grid%dnw, c2a, grid%pb, grid%p, grid%alt, &
1206	grid%msfux,grid%msfuy, grid%msfvx, grid%msfvx_inv, &
1207	grid%msfvy, grid%msftx,grid%msfty, &
1208	grid%rdx, grid%rdy, rk_step, &
1209	ids, ide, jds, jde, kds, kde, &
1210	ims, ime, jms, jme, kms, kme, &
1211	grid%i_start(ij), grid%i_end(ij), &
1212	grid%j_start(ij), grid%j_end(ij), &
1213	k_start , k_end )
1214
1215	#ifdef LMDZ1
1216	WRITE(message, *)' dyn_em: post small_step_prep'
1217	CALL wrf_debug(200, message)
1218	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1219	' u_tend: ', ru_tendf(im2,1,jm2)
1220	CALL wrf_debug(200, message)
1221	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1222	'p sfc: ',p8w(im2,kms,jm2)
1223	CALL wrf_debug(200, message)
1224	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1225	CALL wrf_debug(200, message)
1226	#endif
1227	CALL calc_p_rho( grid%al, grid%p, grid%ph_2, &
1228	grid%alt, grid%t_2, grid%t_save, c2a, pm1, &
1229	grid%mu_2, grid%muts, grid%znu, t0, &
1230	grid%rdnw, grid%dnw, grid%smdiv, &
1231	config_flags%non_hydrostatic, 0, &
1232	ids, ide, jds, jde, kds, kde, &
1233	ims, ime, jms, jme, kms, kme, &
1234	grid%i_start(ij), grid%i_end(ij), &
1235	grid%j_start(ij), grid%j_end(ij), &
1236	k_start , k_end )
1237	#ifdef LMDZ1
1238	WRITE(message, *)' dyn_em: post calc_p_rho'
1239	CALL wrf_debug(200, message)
1240	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1241	' u_tend: ', ru_tendf(im2,1,jm2)
1242	CALL wrf_debug(200, message)
1243	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1244	'p sfc: ',p8w(im2,kms,jm2)
1245	CALL wrf_debug(200, message)
1246	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1247	CALL wrf_debug(200, message)
1248	#endif
1249
1250	IF (config_flags%non_hydrostatic) THEN
1251	CALL calc_coef_w( a,alpha,gamma, &
1252	grid%mut, cqw, &
1253	grid%rdn, grid%rdnw, c2a, &
1254	dts_rk, g, grid%epssm, &
1255	config_flags%top_lid, &
1256	ids, ide, jds, jde, kds, kde, &
1257	ims, ime, jms, jme, kms, kme, &
1258	grid%i_start(ij), grid%i_end(ij), &
1259	grid%j_start(ij), grid%j_end(ij), &
1260	k_start , k_end )
1261	ENDIF
1262	#ifdef LMDZ1
1263	WRITE(message, *)' dyn_em: post calc_coef_w'
1264	CALL wrf_debug(200, message)
1265	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1266	' u_tend: ', ru_tendf(im2,1,jm2)
1267	CALL wrf_debug(200, message)
1268	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1269	'p sfc: ',p8w(im2,kms,jm2)
1270	CALL wrf_debug(200, message)
1271	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1272	CALL wrf_debug(200, message)
1273	WRITE(message,*)' al: ',grid%al(im2,km2,jm2), ' p: ', grid%p(im2,1,jm2), &
1274	' ph: ',grid%ph_2(im2,1,jm2)
1275	CALL wrf_debug(200, message)
1276	WRITE(message,*)' mu: ',grid%mu_2(im2,jm2), 'alt: ', grid%alt(im2,1,jm2), &
1277	' mu: ',grid%mu_2(im2,jm2), ' znu: ', grid%znu(1),' ph 1: ',grid%ph_2(im2,2,jm2)
1278	CALL wrf_debug(200, message)
1279	PRINT *,' c2a: ',c2a(im2,1,jm2), 't_2: ', grid%t_2(im2,1,jm2), &
1280	' t_save: ',grid%t_save(im2,1,jm2),' pm1: ',pm1(im2,1,jm2)
1281	CALL wrf_debug(200, message)
1282	#endif
1283	ENDDO
1284	!$OMP END PARALLEL DO
1285	BENCH_END(small_step_prep_tim)
1286
1287	#ifdef DM_PARALLEL
1288	!-----------------------------------------------------------------------
1289	! Stencils for patch communications (WCS, 29 June 2001)
1290	! Note: the small size of this halo exchange reflects the
1291	! fact that we are carrying the uncoupled variables
1292	! as state variables in the mass coordinate model, as
1293	! opposed to the coupled variables as in the height
1294	! coordinate model.
1295	!
1296	! * * * * *
1297	! * * * * * * * * *
1298	! * + * * + * * * + * *
1299	! * * * * * * * * *
1300	! * * * * *
1301	!
1302	! 3D variables - note staggering! ph_2(Z), u_save(X), v_save(Y)
1303	!
1304	! ph_2 x
1305	! al x
1306	! p x
1307	! t_1 x
1308	! t_save x
1309	! u_save x
1310	! v_save x
1311	!
1312	! the following are 2D (xy) variables
1313	!
1314	! mu_1 x
1315	! mu_2 x
1316	! mudf x
1317	! php x
1318	! alt x
1319	! pb x
1320	!--------------------------------------------------------------
1321	# include "HALO_EM_B.inc"
1322	# include "PERIOD_BDY_EM_B.inc"
1323	#endif
1324
1325	BENCH_START(set_phys_bc2_tim)
1326	!$OMP PARALLEL DO &
1327	!$OMP PRIVATE ( ij )
1328
1329	DO ij = 1 , grid%num_tiles
1330
1331	CALL set_physical_bc3d( grid%ru_tend, 'u', config_flags, &
1332	ids, ide, jds, jde, kds, kde, &
1333	ims, ime, jms, jme, kms, kme, &
1334	ips, ipe, jps, jpe, kps, kpe, &
1335	grid%i_start(ij), grid%i_end(ij), &
1336	grid%j_start(ij), grid%j_end(ij), &
1337	k_start , k_end )
1338
1339	CALL set_physical_bc3d( grid%rv_tend, 'v', config_flags, &
1340	ids, ide, jds, jde, kds, kde, &
1341	ims, ime, jms, jme, kms, kme, &
1342	ips, ipe, jps, jpe, kps, kpe, &
1343	grid%i_start(ij), grid%i_end(ij), &
1344	grid%j_start(ij), grid%j_end(ij), &
1345	k_start , k_end )
1346
1347	CALL set_physical_bc3d( grid%ph_2, 'w', config_flags, &
1348	ids, ide, jds, jde, kds, kde, &
1349	ims, ime, jms, jme, kms, kme, &
1350	ips, ipe, jps, jpe, kps, kpe, &
1351	grid%i_start(ij), grid%i_end(ij), &
1352	grid%j_start(ij), grid%j_end(ij), &
1353	k_start , k_end )
1354
1355	CALL set_physical_bc3d( grid%al, 'p', config_flags, &
1356	ids, ide, jds, jde, kds, kde, &
1357	ims, ime, jms, jme, kms, kme, &
1358	ips, ipe, jps, jpe, kps, kpe, &
1359	grid%i_start(ij), grid%i_end(ij), &
1360	grid%j_start(ij), grid%j_end(ij), &
1361	k_start , k_end )
1362
1363	CALL set_physical_bc3d( grid%p, 'p', config_flags, &
1364	ids, ide, jds, jde, kds, kde, &
1365	ims, ime, jms, jme, kms, kme, &
1366	ips, ipe, jps, jpe, kps, kpe, &
1367	grid%i_start(ij), grid%i_end(ij), &
1368	grid%j_start(ij), grid%j_end(ij), &
1369	k_start , k_end )
1370
1371	CALL set_physical_bc3d( grid%t_1, 'p', config_flags, &
1372	ids, ide, jds, jde, kds, kde, &
1373	ims, ime, jms, jme, kms, kme, &
1374	ips, ipe, jps, jpe, kps, kpe, &
1375	grid%i_start(ij), grid%i_end(ij), &
1376	grid%j_start(ij), grid%j_end(ij), &
1377	k_start , k_end )
1378
1379	CALL set_physical_bc3d( grid%t_save, 't', config_flags, &
1380	ids, ide, jds, jde, kds, kde, &
1381	ims, ime, jms, jme, kms, kme, &
1382	ips, ipe, jps, jpe, kps, kpe, &
1383	grid%i_start(ij), grid%i_end(ij), &
1384	grid%j_start(ij), grid%j_end(ij), &
1385	k_start , k_end )
1386
1387	CALL set_physical_bc2d( grid%mu_1, 't', config_flags, &
1388	ids, ide, jds, jde, &
1389	ims, ime, jms, jme, &
1390	ips, ipe, jps, jpe, &
1391	grid%i_start(ij), grid%i_end(ij), &
1392	grid%j_start(ij), grid%j_end(ij) )
1393
1394	CALL set_physical_bc2d( grid%mu_2, 't', config_flags, &
1395	ids, ide, jds, jde, &
1396	ims, ime, jms, jme, &
1397	ips, ipe, jps, jpe, &
1398	grid%i_start(ij), grid%i_end(ij), &
1399	grid%j_start(ij), grid%j_end(ij) )
1400
1401	CALL set_physical_bc2d( grid%mudf, 't', config_flags, &
1402	ids, ide, jds, jde, &
1403	ims, ime, jms, jme, &
1404	ips, ipe, jps, jpe, &
1405	grid%i_start(ij), grid%i_end(ij), &
1406	grid%j_start(ij), grid%j_end(ij) )
1407
1408	END DO
1409	!$OMP END PARALLEL DO
1410	BENCH_END(set_phys_bc2_tim)
1411	small_steps : DO iteration = 1 , number_of_small_timesteps
1412
1413	! Boundary condition time (or communication time).
1414	#ifdef DM_PARALLEL
1415	# include "PERIOD_BDY_EM_B.inc"
1416	#endif
1417
1418	!$OMP PARALLEL DO &
1419	!$OMP PRIVATE ( ij )
1420
1421	DO ij = 1 , grid%num_tiles
1422
1423	BENCH_START(advance_uv_tim)
1424	#ifdef LMDZ1
1425	WRITE(message, *)' dyn_em: before advance_uv'
1426	CALL wrf_debug(200, message)
1427	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1428	' u_tend: ', ru_tendf(im2,1,jm2)
1429	CALL wrf_debug(200, message)
1430	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1431	'p sfc: ',p8w(im2,kms,jm2)
1432	CALL wrf_debug(200, message)
1433	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1434	CALL wrf_debug(200, message)
1435	WRITE(message,*)' al: ',grid%al(im2,km2,jm2), ' p: ', grid%p(im2,1,jm2), &
1436	' ph: ',grid%ph_2(im2,1,jm2)
1437	CALL wrf_debug(200, message)
1438	WRITE(message,*)' mu: ',grid%mu_2(im2,jm2), 'alt: ', grid%alt(im2,1,jm2), &
1439	' mu: ',grid%mu_2(im2,jm2), ' znu: ', grid%znu(1),' ph 1: ',grid%ph_2(im2,2,jm2)
1440	CALL wrf_debug(200, message)
1441	PRINT *,' c2a: ',c2a(im2,1,jm2), 't_2: ', grid%t_2(im2,1,jm2), &
1442	' t_save: ',grid%t_save(im2,1,jm2),' pm1: ',pm1(im2,1,jm2)
1443	CALL wrf_debug(200, message)
1444	#endif
1445
1446	CALL advance_uv ( grid%u_2, grid%ru_tend, grid%v_2, grid%rv_tend, &
1447	grid%p, grid%pb, &
1448	grid%ph_2, grid%php, grid%alt, grid%al, &
1449	grid%mu_2, &
1450	grid%muu, cqu, grid%muv, cqv, grid%mudf, &
1451	grid%msfux, grid%msfuy, grid%msfvx, &
1452	grid%msfvx_inv, grid%msfvy, &
1453	grid%rdx, grid%rdy, dts_rk, &
1454	grid%cf1, grid%cf2, grid%cf3, grid%fnm, grid%fnp, &
1455	grid%emdiv, &
1456	grid%rdnw, config_flags,grid%spec_zone, &
1457	config_flags%non_hydrostatic, config_flags%top_lid, &
1458	ids, ide, jds, jde, kds, kde, &
1459	ims, ime, jms, jme, kms, kme, &
1460	grid%i_start(ij), grid%i_end(ij), &
1461	grid%j_start(ij), grid%j_end(ij), &
1462	k_start , k_end )
1463	BENCH_END(advance_uv_tim)
1464
1465	END DO
1466	!$OMP END PARALLEL DO
1467	#ifdef LMDZ1
1468	WRITE(message, *)' dyn_em: after advance_uv'
1469	CALL wrf_debug(200, message)
1470	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1471	' u_tend: ', ru_tendf(im2,1,jm2)
1472	CALL wrf_debug(200, message)
1473	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1474	'p sfc: ',p8w(im2,kms,jm2)
1475	CALL wrf_debug(200, message)
1476	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1477	CALL wrf_debug(200, message)
1478	WRITE(message,*)' al: ',grid%al(im2,km2,jm2), ' p: ', grid%p(im2,1,jm2), &
1479	' ph: ',grid%ph_2(im2,1,jm2)
1480	CALL wrf_debug(200, message)
1481	WRITE(message,*)' mu: ',grid%mu_2(im2,jm2), 'alt: ', grid%alt(im2,1,jm2), &
1482	' mu: ',grid%mu_2(im2,jm2), ' znu: ', grid%znu(1),' ph 1: ',grid%ph_2(im2,2,jm2)
1483	CALL wrf_debug(200, message)
1484	PRINT *,' c2a: ',c2a(im2,1,jm2), 't_2: ', grid%t_2(im2,1,jm2), &
1485	' t_save: ',grid%t_save(im2,1,jm2),' pm1: ',pm1(im2,1,jm2)
1486	CALL wrf_debug(200, message)
1487	#endif
1488
1489	!-----------------------------------------------------------
1490	! acoustic integration polar filter for smallstep u, v
1491	!-----------------------------------------------------------
1492
1493	IF (config_flags%polar) THEN
1494
1495	CALL pxft ( grid=grid &
1496	,lineno=__LINE__ &
1497	,flag_uv = 1 &
1498	,flag_rurv = 0 &
1499	,flag_wph = 0 &
1500	,flag_ww = 0 &
1501	,flag_t = 0 &
1502	,flag_mu = 0 &
1503	,flag_mut = 0 &
1504	,flag_moist = 0 &
1505	,flag_chem = 0 &
1506	,flag_tracer = 0 &
1507	,flag_scalar = 0 &
1508	,positive_definite = .FALSE. &
1509	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
1510	,fft_filter_lat = config_flags%fft_filter_lat &
1511	,dclat = dclat &
1512	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
1513	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
1514	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
1515	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
1516	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
1517
1518	END IF
1519
1520	!-----------------------------------------------------------
1521	! end acoustic integration polar filter for smallstep u, v
1522	!-----------------------------------------------------------
1523
1524	!$OMP PARALLEL DO &
1525	!$OMP PRIVATE ( ij )
1526	DO ij = 1 , grid%num_tiles
1527
1528	BENCH_START(spec_bdy_uv_tim)
1529	IF( config_flags%specified .or. config_flags%nested ) THEN
1530	CALL spec_bdyupdate(grid%u_2, grid%ru_tend, dts_rk, &
1531	'u' , config_flags, &
1532	grid%spec_zone, &
1533	ids,ide, jds,jde, kds,kde, & ! domain dims
1534	ims,ime, jms,jme, kms,kme, & ! memory dims
1535	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
1536	grid%i_start(ij), grid%i_end(ij), &
1537	grid%j_start(ij), grid%j_end(ij), &
1538	k_start , k_end )
1539
1540	CALL spec_bdyupdate(grid%v_2, grid%rv_tend, dts_rk, &
1541	'v' , config_flags, &
1542	grid%spec_zone, &
1543	ids,ide, jds,jde, kds,kde, & ! domain dims
1544	ims,ime, jms,jme, kms,kme, & ! memory dims
1545	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
1546	grid%i_start(ij), grid%i_end(ij), &
1547	grid%j_start(ij), grid%j_end(ij), &
1548	k_start , k_end )
1549
1550	ENDIF
1551	BENCH_END(spec_bdy_uv_tim)
1552
1553	END DO
1554	!$OMP END PARALLEL DO
1555
1556	#ifdef DM_PARALLEL
1557	!
1558	! Stencils for patch communications (WCS, 29 June 2001)
1559	!
1560	! * *
1561	! * + * * + * +
1562	! * *
1563	!
1564	! u_2 x
1565	! v_2 x
1566	!
1567	# include "HALO_EM_C.inc"
1568	#endif
1569
1570	!$OMP PARALLEL DO &
1571	!$OMP PRIVATE ( ij )
1572	DO ij = 1 , grid%num_tiles
1573
1574	! advance the mass in the column, theta, and calculate ww
1575
1576	BENCH_START(advance_mu_t_tim)
1577	CALL advance_mu_t( grid%ww, ww1, grid%u_2, grid%u_save, grid%v_2, grid%v_save, &
1578	grid%mu_2, grid%mut, muave, grid%muts, grid%muu, grid%muv, &
1579	grid%mudf, grid%ru_m, grid%rv_m, grid%ww_m, &
1580	grid%t_2, grid%t_save, t_2save, t_tend, &
1581	mu_tend, &
1582	grid%rdx, grid%rdy, dts_rk, grid%epssm, &
1583	grid%dnw, grid%fnm, grid%fnp, grid%rdnw, &
1584	grid%msfux,grid%msfuy, grid%msfvx, grid%msfvx_inv, &
1585	grid%msfvy, grid%msftx,grid%msfty, &
1586	iteration, config_flags, &
1587	ids, ide, jds, jde, kds, kde, &
1588	ims, ime, jms, jme, kms, kme, &
1589	grid%i_start(ij), grid%i_end(ij), &
1590	grid%j_start(ij), grid%j_end(ij), &
1591	k_start , k_end )
1592	BENCH_END(advance_mu_t_tim)
1593	ENDDO
1594	!$OMP END PARALLEL DO
1595	#ifdef LMDZ1
1596	WRITE(message, *)' dyn_em: after advance_mut'
1597	CALL wrf_debug(200, message)
1598	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1599	' u_tend: ', ru_tendf(im2,1,jm2)
1600	CALL wrf_debug(200, message)
1601	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1602	'p sfc: ',p8w(im2,kms,jm2)
1603	CALL wrf_debug(200, message)
1604	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1605	CALL wrf_debug(200, message)
1606	WRITE(message,*)' al: ',grid%al(im2,km2,jm2), ' p: ', grid%p(im2,1,jm2), &
1607	' ph: ',grid%ph_2(im2,1,jm2)
1608	CALL wrf_debug(200, message)
1609	WRITE(message,*)' mu: ',grid%mu_2(im2,jm2), 'alt: ', grid%alt(im2,1,jm2), &
1610	' mu: ',grid%mu_2(im2,jm2), ' znu: ', grid%znu(1),' ph 1: ',grid%ph_2(im2,2,jm2)
1611	CALL wrf_debug(200, message)
1612	PRINT *,' c2a: ',c2a(im2,1,jm2), 't_2: ', grid%t_2(im2,1,jm2), &
1613	' t_save: ',grid%t_save(im2,1,jm2),' pm1: ',pm1(im2,1,jm2)
1614	CALL wrf_debug(200, message)
1615	#endif
1616
1617	!-----------------------------------------------------------
1618	! acoustic integration polar filter for smallstep mu, t
1619	!-----------------------------------------------------------
1620
1621	IF ( (config_flags%polar) ) THEN
1622
1623	CALL pxft ( grid=grid &
1624	,lineno=__LINE__ &
1625	,flag_uv = 0 &
1626	,flag_rurv = 0 &
1627	,flag_wph = 0 &
1628	,flag_ww = 0 &
1629	,flag_t = 1 &
1630	,flag_mu = 1 &
1631	,flag_mut = 0 &
1632	,flag_moist = 0 &
1633	,flag_chem = 0 &
1634	,flag_tracer = 0 &
1635	,flag_scalar = 0 &
1636	,positive_definite = .FALSE. &
1637	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
1638	,fft_filter_lat = config_flags%fft_filter_lat &
1639	,dclat = dclat &
1640	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
1641	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
1642	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
1643	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
1644	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
1645
1646	grid%muts = grid%mut + grid%mu_2 ! reset muts using filtered mu_2
1647
1648	END IF
1649
1650	!-----------------------------------------------------------
1651	! end acoustic integration polar filter for smallstep mu, t
1652	!-----------------------------------------------------------
1653
1654	BENCH_START(spec_bdy_t_tim)
1655
1656	!$OMP PARALLEL DO &
1657	!$OMP PRIVATE ( ij )
1658	DO ij = 1 , grid%num_tiles
1659
1660	IF( config_flags%specified .or. config_flags%nested ) THEN
1661
1662	CALL spec_bdyupdate(grid%t_2, t_tend, dts_rk, &
1663	't' , config_flags, &
1664	grid%spec_zone, &
1665	ids,ide, jds,jde, kds,kde, &
1666	ims,ime, jms,jme, kms,kme, &
1667	ips,ipe, jps,jpe, kps,kpe, &
1668	grid%i_start(ij), grid%i_end(ij),&
1669	grid%j_start(ij), grid%j_end(ij),&
1670	k_start , k_end )
1671
1672	CALL spec_bdyupdate(grid%mu_2, mu_tend, dts_rk, &
1673	'm' , config_flags, &
1674	grid%spec_zone, &
1675	ids,ide, jds,jde, 1 ,1 , &
1676	ims,ime, jms,jme, 1 ,1 , &
1677	ips,ipe, jps,jpe, 1 ,1 , &
1678	grid%i_start(ij), grid%i_end(ij),&
1679	grid%j_start(ij), grid%j_end(ij),&
1680	1 , 1 )
1681
1682	CALL spec_bdyupdate(grid%muts, mu_tend, dts_rk, &
1683	'm' , config_flags, &
1684	grid%spec_zone, &
1685	ids,ide, jds,jde, 1 ,1 , & ! domain dims
1686	ims,ime, jms,jme, 1 ,1 , & ! memory dims
1687	ips,ipe, jps,jpe, 1 ,1 , & ! patch dims
1688	grid%i_start(ij), grid%i_end(ij), &
1689	grid%j_start(ij), grid%j_end(ij), &
1690	1 , 1 )
1691	ENDIF
1692	BENCH_END(spec_bdy_t_tim)
1693
1694	! small (acoustic) step for the vertical momentum,
1695	! density and coupled potential temperature.
1696
1697
1698	BENCH_START(advance_w_tim)
1699	IF ( config_flags%non_hydrostatic ) THEN
1700	CALL advance_w( grid%w_2, rw_tend, grid%ww, w_save, &
1701	grid%u_2, grid%v_2, &
1702	grid%mu_2, grid%mut, muave, grid%muts, &
1703	t_2save, grid%t_2, grid%t_save, &
1704	grid%ph_2, ph_save, grid%phb, ph_tend, &
1705	grid%ht, c2a, cqw, grid%alt, grid%alb, &
1706	a, alpha, gamma, &
1707	grid%rdx, grid%rdy, dts_rk, t0, grid%epssm, &
1708	grid%dnw, grid%fnm, grid%fnp, grid%rdnw, &
1709	grid%rdn, grid%cf1, grid%cf2, grid%cf3, &
1710	grid%msftx, grid%msfty, &
1711	config_flags, config_flags%top_lid, &
1712	ids,ide, jds,jde, kds,kde, &
1713	ims,ime, jms,jme, kms,kme, &
1714	grid%i_start(ij), grid%i_end(ij), &
1715	grid%j_start(ij), grid%j_end(ij), &
1716	k_start , k_end )
1717	ENDIF
1718	BENCH_END(advance_w_tim)
1719
1720	ENDDO
1721	!$OMP END PARALLEL DO
1722	#ifdef LMDZ1
1723	WRITE(message, *)' dyn_em: after advance_w'
1724	CALL wrf_debug(200, message)
1725	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1726	' u_tend: ', ru_tendf(im2,1,jm2)
1727	CALL wrf_debug(200, message)
1728	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1729	'p sfc: ',p8w(im2,kms,jm2)
1730	CALL wrf_debug(200, message)
1731	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1732	CALL wrf_debug(200, message)
1733	WRITE(message,*)' al: ',grid%al(im2,km2,jm2), ' p: ', grid%p(im2,1,jm2), &
1734	' ph: ',grid%ph_2(im2,1,jm2)
1735	CALL wrf_debug(200, message)
1736	WRITE(message,*)' mu: ',grid%mu_2(im2,jm2), 'alt: ', grid%alt(im2,1,jm2), &
1737	' mu: ',grid%mu_2(im2,jm2), ' znu: ', grid%znu(1),' ph 1: ',grid%ph_2(im2,2,jm2)
1738	CALL wrf_debug(200, message)
1739	PRINT *,' c2a: ',c2a(im2,1,jm2), 't_2: ', grid%t_2(im2,1,jm2), &
1740	' t_save: ',grid%t_save(im2,1,jm2),' pm1: ',pm1(im2,1,jm2)
1741	CALL wrf_debug(200, message)
1742	#endif
1743
1744	!-----------------------------------------------------------
1745	! acoustic integration polar filter for smallstep w, geopotential
1746	!-----------------------------------------------------------
1747
1748	IF ( (config_flags%polar) .AND. (config_flags%non_hydrostatic) ) THEN
1749
1750	CALL pxft ( grid=grid &
1751	,lineno=__LINE__ &
1752	,flag_uv = 0 &
1753	,flag_rurv = 0 &
1754	,flag_wph = 1 &
1755	,flag_ww = 0 &
1756	,flag_t = 0 &
1757	,flag_mu = 0 &
1758	,flag_mut = 0 &
1759	,flag_moist = 0 &
1760	,flag_chem = 0 &
1761	,flag_tracer = 0 &
1762	,flag_scalar = 0 &
1763	,positive_definite = .FALSE. &
1764	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
1765	,fft_filter_lat = config_flags%fft_filter_lat &
1766	,dclat = dclat &
1767	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
1768	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
1769	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
1770	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
1771	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
1772
1773	END IF
1774
1775	!-----------------------------------------------------------
1776	! end acoustic integration polar filter for smallstep w, geopotential
1777	!-----------------------------------------------------------
1778
1779	!$OMP PARALLEL DO &
1780	!$OMP PRIVATE ( ij )
1781	DO ij = 1 , grid%num_tiles
1782
1783	BENCH_START(sumflux_tim)
1784	CALL sumflux ( grid%u_2, grid%v_2, grid%ww, &
1785	grid%u_save, grid%v_save, ww1, &
1786	grid%muu, grid%muv, &
1787	grid%ru_m, grid%rv_m, grid%ww_m, grid%epssm, &
1788	grid%msfux, grid% msfuy, grid%msfvx, &
1789	grid%msfvx_inv, grid%msfvy, &
1790	iteration, number_of_small_timesteps, &
1791	ids, ide, jds, jde, kds, kde, &
1792	ims, ime, jms, jme, kms, kme, &
1793	grid%i_start(ij), grid%i_end(ij), &
1794	grid%j_start(ij), grid%j_end(ij), &
1795	k_start , k_end )
1796	BENCH_END(sumflux_tim)
1797	#ifdef LMDZ1
1798	WRITE(message, *)' dyn_em: after sumflux'
1799	CALL wrf_debug(200, message)
1800	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1801	' u_tend: ', ru_tendf(im2,1,jm2)
1802	CALL wrf_debug(200, message)
1803	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1804	'p sfc: ',p8w(im2,kms,jm2)
1805	CALL wrf_debug(200, message)
1806	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1807	CALL wrf_debug(200, message)
1808	WRITE(message,*)' al: ',grid%al(im2,km2,jm2), ' p: ', grid%p(im2,1,jm2), &
1809	' ph: ',grid%ph_2(im2,1,jm2)
1810	CALL wrf_debug(200, message)
1811	WRITE(message,*)' mu: ',grid%mu_2(im2,jm2), 'alt: ', grid%alt(im2,1,jm2), &
1812	' mu: ',grid%mu_2(im2,jm2), ' znu: ', grid%znu(1),' ph 1: ',grid%ph_2(im2,2,jm2)
1813	CALL wrf_debug(200, message)
1814	PRINT *,' c2a: ',c2a(im2,1,jm2), 't_2: ', grid%t_2(im2,1,jm2), &
1815	' t_save: ',grid%t_save(im2,1,jm2),' pm1: ',pm1(im2,1,jm2)
1816	CALL wrf_debug(200, message)
1817	#endif
1818
1819	IF( config_flags%specified .or. config_flags%nested ) THEN
1820
1821	BENCH_START(spec_bdynhyd_tim)
1822	IF (config_flags%non_hydrostatic) THEN
1823	CALL spec_bdyupdate_ph( ph_save, grid%ph_2, ph_tend, &
1824	mu_tend, grid%muts, dts_rk, &
1825	'h' , config_flags, &
1826	grid%spec_zone, &
1827	ids,ide, jds,jde, kds,kde, &
1828	ims,ime, jms,jme, kms,kme, &
1829	ips,ipe, jps,jpe, kps,kpe, &
1830	grid%i_start(ij), grid%i_end(ij),&
1831	grid%j_start(ij), grid%j_end(ij),&
1832	k_start , k_end )
1833	#ifdef LMDZ1
1834	WRITE(message, *)' dyn_em: after spec_bdynhyd_ph'
1835	CALL wrf_debug(200, message)
1836	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1837	' u_tend: ', ru_tendf(im2,1,jm2)
1838	CALL wrf_debug(200, message)
1839	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1840	'p sfc: ',p8w(im2,kms,jm2)
1841	CALL wrf_debug(200, message)
1842	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1843	CALL wrf_debug(200, message)
1844	WRITE(message,*)' al: ',grid%al(im2,km2,jm2), ' p: ', grid%p(im2,1,jm2), &
1845	' ph: ',grid%ph_2(im2,1,jm2)
1846	CALL wrf_debug(200, message)
1847	WRITE(message,*)' mu: ',grid%mu_2(im2,jm2), 'alt: ', grid%alt(im2,1,jm2), &
1848	' mu: ',grid%mu_2(im2,jm2), ' znu: ', grid%znu(1),' ph 1: ',grid%ph_2(im2,2,jm2)
1849	CALL wrf_debug(200, message)
1850	PRINT *,' c2a: ',c2a(im2,1,jm2), 't_2: ', grid%t_2(im2,1,jm2), &
1851	' t_save: ',grid%t_save(im2,1,jm2),' pm1: ',pm1(im2,1,jm2)
1852	CALL wrf_debug(200, message)
1853	#endif
1854	IF( config_flags%specified ) THEN
1855	CALL zero_grad_bdy ( grid%w_2, &
1856	'w' , config_flags, &
1857	grid%spec_zone, &
1858	ids,ide, jds,jde, kds,kde, &
1859	ims,ime, jms,jme, kms,kme, &
1860	ips,ipe, jps,jpe, kps,kpe, &
1861	grid%i_start(ij), grid%i_end(ij), &
1862	grid%j_start(ij), grid%j_end(ij), &
1863	k_start , k_end )
1864	ELSE
1865	CALL spec_bdyupdate ( grid%w_2, rw_tend, dts_rk, &
1866	'h' , config_flags, &
1867	grid%spec_zone, &
1868	ids,ide, jds,jde, kds,kde, &
1869	ims,ime, jms,jme, kms,kme, &
1870	ips,ipe, jps,jpe, kps,kpe, &
1871	grid%i_start(ij), grid%i_end(ij),&
1872	grid%j_start(ij), grid%j_end(ij),&
1873	k_start , k_end )
1874	ENDIF
1875	ENDIF
1876	BENCH_END(spec_bdynhyd_tim)
1877	ENDIF
1878	#ifdef LMDZ1
1879	WRITE(message, *)' dyn_em: after spec_bdynhyd_tim'
1880	CALL wrf_debug(200, message)
1881	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1882	' u_tend: ', ru_tendf(im2,1,jm2)
1883	CALL wrf_debug(200, message)
1884	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1885	'p sfc: ',p8w(im2,kms,jm2)
1886	CALL wrf_debug(200, message)
1887	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1888	CALL wrf_debug(200, message)
1889	WRITE(message,*)' al: ',grid%al(im2,km2,jm2), ' p: ', grid%p(im2,1,jm2), &
1890	' ph: ',grid%ph_2(im2,1,jm2)
1891	CALL wrf_debug(200, message)
1892	WRITE(message,*)' mu: ',grid%mu_2(im2,jm2), 'alt: ', grid%alt(im2,1,jm2), &
1893	' mu: ',grid%mu_2(im2,jm2), ' znu: ', grid%znu(1),' ph 1: ',grid%ph_2(im2,2,jm2)
1894	CALL wrf_debug(200, message)
1895	PRINT *,' c2a: ',c2a(im2,1,jm2), 't_2: ', grid%t_2(im2,1,jm2), &
1896	' t_save: ',grid%t_save(im2,1,jm2),' pm1: ',pm1(im2,1,jm2)
1897	CALL wrf_debug(200, message)
1898	#endif
1899
1900	BENCH_START(cald_p_rho_tim)
1901	CALL calc_p_rho( grid%al, grid%p, grid%ph_2, &
1902	grid%alt, grid%t_2, grid%t_save, c2a, pm1, &
1903	grid%mu_2, grid%muts, grid%znu, t0, &
1904	grid%rdnw, grid%dnw, grid%smdiv, &
1905	config_flags%non_hydrostatic, iteration, &
1906	ids, ide, jds, jde, kds, kde, &
1907	ims, ime, jms, jme, kms, kme, &
1908	grid%i_start(ij), grid%i_end(ij), &
1909	grid%j_start(ij), grid%j_end(ij), &
1910	k_start , k_end )
1911	BENCH_END(cald_p_rho_tim)
1912	#ifdef LMDZ1
1913	WRITE(message, *)' dyn_em: after calc_p_rho'
1914	CALL wrf_debug(200, message)
1915	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1916	' u_tend: ', ru_tendf(im2,1,jm2)
1917	CALL wrf_debug(200, message)
1918	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1919	'p sfc: ',p8w(im2,kms,jm2)
1920	CALL wrf_debug(200, message)
1921	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1922	CALL wrf_debug(200, message)
1923	#endif
1924
1925	ENDDO
1926	!$OMP END PARALLEL DO
1927	#ifdef LMDZ1
1928	WRITE(message, *)' dyn_em: after geopotential'
1929	CALL wrf_debug(200, message)
1930	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
1931	' u_tend: ', ru_tendf(im2,1,jm2)
1932	CALL wrf_debug(200, message)
1933	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
1934	'p sfc: ',p8w(im2,kms,jm2)
1935	CALL wrf_debug(200, message)
1936	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
1937	CALL wrf_debug(200, message)
1938	#endif
1939
1940	#ifdef DM_PARALLEL
1941	!
1942	! Stencils for patch communications (WCS, 29 June 2001)
1943	!
1944	! * *
1945	! * + * * + * +
1946	! * *
1947	!
1948	! ph_2 x
1949	! al x
1950	! p x
1951	!
1952	! 2D variables (x,y)
1953	!
1954	! mu_2 x
1955	! muts x
1956	! mudf x
1957
1958	# include "HALO_EM_C2.inc"
1959	# include "PERIOD_BDY_EM_B3.inc"
1960	#endif
1961
1962	BENCH_START(phys_bc_tim)
1963	!$OMP PARALLEL DO &
1964	!$OMP PRIVATE ( ij )
1965	DO ij = 1 , grid%num_tiles
1966
1967	! boundary condition set for next small timestep
1968
1969	CALL set_physical_bc3d( grid%ph_2, 'w', config_flags, &
1970	ids, ide, jds, jde, kds, kde, &
1971	ims, ime, jms, jme, kms, kme, &
1972	ips, ipe, jps, jpe, kps, kpe, &
1973	grid%i_start(ij), grid%i_end(ij), &
1974	grid%j_start(ij), grid%j_end(ij), &
1975	k_start , k_end )
1976
1977	CALL set_physical_bc3d( grid%al, 'p', config_flags, &
1978	ids, ide, jds, jde, kds, kde, &
1979	ims, ime, jms, jme, kms, kme, &
1980	ips, ipe, jps, jpe, kps, kpe, &
1981	grid%i_start(ij), grid%i_end(ij), &
1982	grid%j_start(ij), grid%j_end(ij), &
1983	k_start , k_end )
1984
1985	CALL set_physical_bc3d( grid%p, 'p', config_flags, &
1986	ids, ide, jds, jde, kds, kde, &
1987	ims, ime, jms, jme, kms, kme, &
1988	ips, ipe, jps, jpe, kps, kpe, &
1989	grid%i_start(ij), grid%i_end(ij), &
1990	grid%j_start(ij), grid%j_end(ij), &
1991	k_start , k_end )
1992
1993	CALL set_physical_bc2d( grid%muts, 't', config_flags, &
1994	ids, ide, jds, jde, &
1995	ims, ime, jms, jme, &
1996	ips, ipe, jps, jpe, &
1997	grid%i_start(ij), grid%i_end(ij), &
1998	grid%j_start(ij), grid%j_end(ij) )
1999
2000	CALL set_physical_bc2d( grid%mu_2, 't', config_flags, &
2001	ids, ide, jds, jde, &
2002	ims, ime, jms, jme, &
2003	ips, ipe, jps, jpe, &
2004	grid%i_start(ij), grid%i_end(ij), &
2005	grid%j_start(ij), grid%j_end(ij) )
2006
2007	CALL set_physical_bc2d( grid%mudf, 't', config_flags, &
2008	ids, ide, jds, jde, &
2009	ims, ime, jms, jme, &
2010	ips, ipe, jps, jpe, &
2011	grid%i_start(ij), grid%i_end(ij), &
2012	grid%j_start(ij), grid%j_end(ij) )
2013
2014	END DO
2015	!$OMP END PARALLEL DO
2016	BENCH_END(phys_bc_tim)
2017
2018	END DO small_steps
2019
2020	!$OMP PARALLEL DO &
2021	!$OMP PRIVATE ( ij )
2022	DO ij = 1 , grid%num_tiles
2023
2024	CALL wrf_debug ( 200 , ' call rk_small_finish' )
2025
2026	! change time-perturbation variables back to
2027	! full perturbation variables.
2028	! first get updated mu at u and v points
2029
2030	BENCH_START(calc_mu_uv_tim)
2031	CALL calc_mu_uv_1 ( config_flags, &
2032	grid%muts, muus, muvs, &
2033	ids, ide, jds, jde, kds, kde, &
2034	ims, ime, jms, jme, kms, kme, &
2035	grid%i_start(ij), grid%i_end(ij), &
2036	grid%j_start(ij), grid%j_end(ij), &
2037	k_start , k_end )
2038	BENCH_END(calc_mu_uv_tim)
2039	BENCH_START(small_step_finish_tim)
2040	CALL small_step_finish( grid%u_2, grid%u_1, grid%v_2, grid%v_1, grid%w_2, grid%w_1, &
2041	grid%t_2, grid%t_1, grid%ph_2, grid%ph_1, grid%ww, ww1, &
2042	grid%mu_2, grid%mu_1, &
2043	grid%mut, grid%muts, grid%muu, muus, grid%muv, muvs, &
2044	grid%u_save, grid%v_save, w_save, &
2045	grid%t_save, ph_save, mu_save, &
2046	grid%msfux,grid%msfuy, grid%msfvx,grid%msfvy, grid%msftx,grid%msfty, &
2047	grid%h_diabatic, &
2048	number_of_small_timesteps,dts_rk, &
2049	rk_step, rk_order, &
2050	ids, ide, jds, jde, kds, kde, &
2051	ims, ime, jms, jme, kms, kme, &
2052	grid%i_start(ij), grid%i_end(ij), &
2053	grid%j_start(ij), grid%j_end(ij), &
2054	k_start , k_end )
2055	! call to set ru_m, rv_m and ww_m b.c's for PD advection
2056
2057	IF (rk_step == rk_order) THEN
2058
2059	CALL set_physical_bc3d( grid%ru_m, 'u', config_flags, &
2060	ids, ide, jds, jde, kds, kde, &
2061	ims, ime, jms, jme, kms, kme, &
2062	ips, ipe, jps, jpe, kps, kpe, &
2063	grid%i_start(ij), grid%i_end(ij), &
2064	grid%j_start(ij), grid%j_end(ij), &
2065	k_start , k_end )
2066
2067	CALL set_physical_bc3d( grid%rv_m, 'v', config_flags, &
2068	ids, ide, jds, jde, kds, kde, &
2069	ims, ime, jms, jme, kms, kme, &
2070	ips, ipe, jps, jpe, kps, kpe, &
2071	grid%i_start(ij), grid%i_end(ij), &
2072	grid%j_start(ij), grid%j_end(ij), &
2073	k_start , k_end )
2074
2075	CALL set_physical_bc3d( grid%ww_m, 'w', config_flags, &
2076	ids, ide, jds, jde, kds, kde, &
2077	ims, ime, jms, jme, kms, kme, &
2078	ips, ipe, jps, jpe, kps, kpe, &
2079	grid%i_start(ij), grid%i_end(ij), &
2080	grid%j_start(ij), grid%j_end(ij), &
2081	k_start , k_end )
2082
2083	CALL set_physical_bc2d( grid%mut, 't', config_flags, &
2084	ids, ide, jds, jde, &
2085	ims, ime, jms, jme, &
2086	ips, ipe, jps, jpe, &
2087	grid%i_start(ij), grid%i_end(ij), &
2088	grid%j_start(ij), grid%j_end(ij) )
2089
2090	CALL set_physical_bc2d( grid%muts, 't', config_flags, &
2091	ids, ide, jds, jde, &
2092	ims, ime, jms, jme, &
2093	ips, ipe, jps, jpe, &
2094	grid%i_start(ij), grid%i_end(ij), &
2095	grid%j_start(ij), grid%j_end(ij) )
2096
2097	END IF
2098
2099	BENCH_END(small_step_finish_tim)
2100
2101	END DO
2102	!$OMP END PARALLEL DO
2103	#ifdef LMDZ1
2104	WRITE(message, *)' dyn_em: after rk_small_finish'
2105	CALL wrf_debug(200, message)
2106	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
2107	' u_tend: ', ru_tendf(im2,1,jm2)
2108	CALL wrf_debug(200, message)
2109	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
2110	'p sfc: ',p8w(im2,kms,jm2)
2111	CALL wrf_debug(200, message)
2112	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
2113	CALL wrf_debug(200, message)
2114	#endif
2115
2116	!-----------------------------------------------------------
2117	! polar filter for full dynamics variables and time-averaged mass fluxes
2118	!-----------------------------------------------------------
2119
2120	IF (config_flags%polar) THEN
2121
2122	CALL pxft ( grid=grid &
2123	,lineno=__LINE__ &
2124	,flag_uv = 1 &
2125	,flag_rurv = 1 &
2126	,flag_wph = 1 &
2127	,flag_ww = 1 &
2128	,flag_t = 1 &
2129	,flag_mu = 1 &
2130	,flag_mut = 1 &
2131	,flag_moist = 0 &
2132	,flag_chem = 0 &
2133	,flag_tracer = 0 &
2134	,flag_scalar = 0 &
2135	,positive_definite = .FALSE. &
2136	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
2137	,fft_filter_lat = config_flags%fft_filter_lat &
2138	,dclat = dclat &
2139	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
2140	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
2141	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
2142	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
2143	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
2144
2145	END IF
2146
2147	!-----------------------------------------------------------
2148	! end polar filter for full dynamics variables and time-averaged mass fluxes
2149	!-----------------------------------------------------------
2150
2151	!-----------------------------------------------------------------------
2152	! add in physics tendency first if positive definite advection is used.
2153	! pd advection applies advective flux limiter on last runge-kutta step
2154	!-----------------------------------------------------------------------
2155	! first moisture
2156
2157	IF ((config_flags%moist_adv_opt /= ORIGINAL) .and. (rk_step == rk_order)) THEN
2158
2159	!$OMP PARALLEL DO &
2160	!$OMP PRIVATE ( ij )
2161	DO ij = 1 , grid%num_tiles
2162	CALL wrf_debug ( 200 , ' call rk_update_scalar_pd' )
2163	DO im = PARAM_FIRST_SCALAR, num_3d_m
2164	CALL rk_update_scalar_pd( im, im, &
2165	moist_old(ims,kms,jms,im), &
2166	moist_tend(ims,kms,jms,im), &
2167	grid%mu_1, grid%mu_1, grid%mub, &
2168	rk_step, dt_rk, grid%spec_zone, &
2169	config_flags, &
2170	ids, ide, jds, jde, kds, kde, &
2171	ims, ime, jms, jme, kms, kme, &
2172	grid%i_start(ij), grid%i_end(ij), &
2173	grid%j_start(ij), grid%j_end(ij), &
2174	k_start , k_end )
2175	ENDDO
2176	END DO
2177	!$OMP END PARALLEL DO
2178
2179	!---------------------- positive definite bc call
2180	#ifdef DM_PARALLEL
2181	IF (config_flags%moist_adv_opt /= ORIGINAL) THEN
2182	IF ( config_flags%h_sca_adv_order <= 4 ) THEN
2183	# include "HALO_EM_MOIST_OLD_E_5.inc"
2184	ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN
2185	# include "HALO_EM_MOIST_OLD_E_7.inc"
2186	ELSE
2187	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
2188	CALL wrf_error_fatal(TRIM(wrf_err_message))
2189	ENDIF
2190	ENDIF
2191	#endif
2192
2193	#ifdef DM_PARALLEL
2194	# include "PERIOD_BDY_EM_MOIST_OLD.inc"
2195	#endif
2196
2197	!$OMP PARALLEL DO &
2198	!$OMP PRIVATE ( ij )
2199	DO ij = 1 , grid%num_tiles
2200	IF (num_3d_m >= PARAM_FIRST_SCALAR) THEN
2201	DO im = PARAM_FIRST_SCALAR , num_3d_m
2202	CALL set_physical_bc3d( moist_old(ims,kms,jms,im), 'p', config_flags, &
2203	ids, ide, jds, jde, kds, kde, &
2204	ims, ime, jms, jme, kms, kme, &
2205	ips, ipe, jps, jpe, kps, kpe, &
2206	grid%i_start(ij), grid%i_end(ij), &
2207	grid%j_start(ij), grid%j_end(ij), &
2208	k_start , k_end )
2209	END DO
2210	ENDIF
2211	END DO
2212	!$OMP END PARALLEL DO
2213
2214	END IF ! end if for moist_adv_opt
2215
2216	! scalars
2217
2218	IF ((config_flags%scalar_adv_opt /= ORIGINAL) .and. (rk_step == rk_order)) THEN
2219
2220	!$OMP PARALLEL DO &
2221	!$OMP PRIVATE ( ij )
2222	DO ij = 1 , grid%num_tiles
2223	CALL wrf_debug ( 200 , ' call rk_update_scalar_pd' )
2224	DO im = PARAM_FIRST_SCALAR, num_3d_s
2225	CALL rk_update_scalar_pd( im, im, &
2226	scalar_old(ims,kms,jms,im), &
2227	scalar_tend(ims,kms,jms,im), &
2228	grid%mu_1, grid%mu_1, grid%mub, &
2229	rk_step, dt_rk, grid%spec_zone, &
2230	config_flags, &
2231	ids, ide, jds, jde, kds, kde, &
2232	ims, ime, jms, jme, kms, kme, &
2233	grid%i_start(ij), grid%i_end(ij), &
2234	grid%j_start(ij), grid%j_end(ij), &
2235	k_start , k_end )
2236	ENDDO
2237	ENDDO
2238	!$OMP END PARALLEL DO
2239
2240	!---------------------- positive definite bc call
2241	#ifdef DM_PARALLEL
2242	IF (config_flags%scalar_adv_opt /= ORIGINAL) THEN
2243	#ifndef RSL
2244	IF ( config_flags%h_sca_adv_order <= 4 ) THEN
2245	# include "HALO_EM_SCALAR_OLD_E_5.inc"
2246	ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN
2247	# include "HALO_EM_SCALAR_OLD_E_7.inc"
2248	ELSE
2249	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
2250	CALL wrf_error_fatal(TRIM(wrf_err_message))
2251	ENDIF
2252	#else
2253	WRITE(wrf_err_message,*)'cannot use pd scheme with RSL - use RSL-LITE'
2254	CALL wrf_error_fatal(TRIM(wrf_err_message))
2255	#endif
2256	endif
2257	#endif
2258
2259	#ifdef DM_PARALLEL
2260	# include "PERIOD_BDY_EM_SCALAR_OLD.inc"
2261	#endif
2262	!$OMP PARALLEL DO &
2263	!$OMP PRIVATE ( ij )
2264
2265	DO ij = 1 , grid%num_tiles
2266	IF (num_3d_s >= PARAM_FIRST_SCALAR) THEN
2267	DO im = PARAM_FIRST_SCALAR , num_3d_s
2268	CALL set_physical_bc3d( scalar_old(ims,kms,jms,im), 'p', config_flags, &
2269	ids, ide, jds, jde, kds, kde, &
2270	ims, ime, jms, jme, kms, kme, &
2271	ips, ipe, jps, jpe, kps, kpe, &
2272	grid%i_start(ij), grid%i_end(ij), &
2273	grid%j_start(ij), grid%j_end(ij), &
2274	k_start , k_end )
2275	END DO
2276	ENDIF
2277	END DO
2278	!$OMP END PARALLEL DO
2279
2280	END IF ! end if for scalar_adv_opt
2281
2282	! chem
2283
2284	IF ((config_flags%chem_adv_opt /= ORIGINAL) .and. (rk_step == rk_order)) THEN
2285
2286	!$OMP PARALLEL DO &
2287	!$OMP PRIVATE ( ij )
2288	DO ij = 1 , grid%num_tiles
2289	CALL wrf_debug ( 200 , ' call rk_update_scalar_pd' )
2290	DO im = PARAM_FIRST_SCALAR, num_3d_c
2291	CALL rk_update_scalar_pd( im, im, &
2292	chem_old(ims,kms,jms,im), &
2293	chem_tend(ims,kms,jms,im), &
2294	grid%mu_1, grid%mu_1, grid%mub, &
2295	rk_step, dt_rk, grid%spec_zone, &
2296	config_flags, &
2297	ids, ide, jds, jde, kds, kde, &
2298	ims, ime, jms, jme, kms, kme, &
2299	grid%i_start(ij), grid%i_end(ij), &
2300	grid%j_start(ij), grid%j_end(ij), &
2301	k_start , k_end )
2302	ENDDO
2303	END DO
2304	!$OMP END PARALLEL DO
2305
2306	!---------------------- positive definite bc call
2307	#ifdef DM_PARALLEL
2308	IF (config_flags%chem_adv_opt /= ORIGINAL) THEN
2309	IF ( config_flags%h_sca_adv_order <= 4 ) THEN
2310	# include "HALO_EM_CHEM_OLD_E_5.inc"
2311	ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN
2312	# include "HALO_EM_CHEM_OLD_E_7.inc"
2313	ELSE
2314	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
2315	CALL wrf_error_fatal(TRIM(wrf_err_message))
2316	ENDIF
2317	ENDIF
2318	#endif
2319
2320	#ifdef DM_PARALLEL
2321	# include "PERIOD_BDY_EM_CHEM_OLD.inc"
2322	#endif
2323
2324	!$OMP PARALLEL DO &
2325	!$OMP PRIVATE ( ij )
2326	DO ij = 1 , grid%num_tiles
2327	IF (num_3d_c >= PARAM_FIRST_SCALAR) THEN
2328	DO im = PARAM_FIRST_SCALAR , num_3d_c
2329	CALL set_physical_bc3d( chem_old(ims,kms,jms,im), 'p', config_flags, &
2330	ids, ide, jds, jde, kds, kde, &
2331	ims, ime, jms, jme, kms, kme, &
2332	ips, ipe, jps, jpe, kps, kpe, &
2333	grid%i_start(ij), grid%i_end(ij), &
2334	grid%j_start(ij), grid%j_end(ij), &
2335	k_start , k_end )
2336	END DO
2337	ENDIF
2338	END DO
2339	!$OMP END PARALLEL DO
2340
2341	ENDIF ! end if for chem_adv_opt
2342
2343	! tracer
2344
2345	IF ((config_flags%tracer_adv_opt /= ORIGINAL) .and. (rk_step == rk_order)) THEN
2346
2347	!$OMP PARALLEL DO &
2348	!$OMP PRIVATE ( ij )
2349	DO ij = 1 , grid%num_tiles
2350	CALL wrf_debug ( 200 , ' call rk_update_scalar_pd' )
2351	DO im = PARAM_FIRST_SCALAR, num_tracer
2352	CALL rk_update_scalar_pd( im, im, &
2353	tracer_old(ims,kms,jms,im), &
2354	tracer_tend(ims,kms,jms,im), &
2355	grid%mu_1, grid%mu_1, grid%mub, &
2356	rk_step, dt_rk, grid%spec_zone, &
2357	config_flags, &
2358	ids, ide, jds, jde, kds, kde, &
2359	ims, ime, jms, jme, kms, kme, &
2360	grid%i_start(ij), grid%i_end(ij), &
2361	grid%j_start(ij), grid%j_end(ij), &
2362	k_start , k_end )
2363	ENDDO
2364	END DO
2365	!$OMP END PARALLEL DO
2366
2367	!---------------------- positive definite bc call
2368	#ifdef DM_PARALLEL
2369	IF (config_flags%tracer_adv_opt /= ORIGINAL) THEN
2370	IF ( config_flags%h_sca_adv_order <= 4 ) THEN
2371	# include "HALO_EM_TRACER_OLD_E_5.inc"
2372	ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN
2373	# include "HALO_EM_TRACER_OLD_E_7.inc"
2374	ELSE
2375	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
2376	CALL wrf_error_fatal(TRIM(wrf_err_message))
2377	ENDIF
2378	ENDIF
2379	#endif
2380
2381	#ifdef DM_PARALLEL
2382	# include "PERIOD_BDY_EM_TRACER_OLD.inc"
2383	#endif
2384
2385	!$OMP PARALLEL DO &
2386	!$OMP PRIVATE ( ij )
2387	DO ij = 1 , grid%num_tiles
2388	IF (num_tracer >= PARAM_FIRST_SCALAR) THEN
2389	DO im = PARAM_FIRST_SCALAR , num_tracer
2390	CALL set_physical_bc3d( tracer_old(ims,kms,jms,im), 'p', config_flags, &
2391	ids, ide, jds, jde, kds, kde, &
2392	ims, ime, jms, jme, kms, kme, &
2393	ips, ipe, jps, jpe, kps, kpe, &
2394	grid%i_start(ij), grid%i_end(ij), &
2395	grid%j_start(ij), grid%j_end(ij), &
2396	k_start , k_end )
2397	END DO
2398	ENDIF
2399	END DO
2400	!$OMP END PARALLEL DO
2401
2402	ENDIF ! end if for tracer_adv_opt
2403
2404	! tke
2405
2406	IF ((config_flags%tke_adv_opt /= ORIGINAL) .and. (rk_step == rk_order) &
2407	.and. (config_flags%km_opt .eq. 2) ) THEN
2408
2409	!$OMP PARALLEL DO &
2410	!$OMP PRIVATE ( ij )
2411	DO ij = 1 , grid%num_tiles
2412	CALL wrf_debug ( 200 , ' call rk_update_scalar_pd' )
2413	CALL rk_update_scalar_pd( 1, 1, &
2414	grid%tke_1, &
2415	tke_tend(ims,kms,jms), &
2416	grid%mu_1, grid%mu_1, grid%mub, &
2417	rk_step, dt_rk, grid%spec_zone, &
2418	config_flags, &
2419	ids, ide, jds, jde, kds, kde, &
2420	ims, ime, jms, jme, kms, kme, &
2421	grid%i_start(ij), grid%i_end(ij), &
2422	grid%j_start(ij), grid%j_end(ij), &
2423	k_start , k_end )
2424	ENDDO
2425	!$OMP END PARALLEL DO
2426
2427	!---------------------- positive definite bc call
2428	#ifdef DM_PARALLEL
2429	IF (config_flags%tke_adv_opt /= ORIGINAL) THEN
2430	IF ( config_flags%h_sca_adv_order <= 4 ) THEN
2431	# include "HALO_EM_TKE_OLD_E_5.inc"
2432	ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN
2433	# include "HALO_EM_TKE_OLD_E_7.inc"
2434	ELSE
2435	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
2436	CALL wrf_error_fatal(TRIM(wrf_err_message))
2437	ENDIF
2438	ENDIF
2439	#endif
2440
2441	#ifdef DM_PARALLEL
2442	# include "PERIOD_BDY_EM_TKE_OLD.inc"
2443	#endif
2444
2445	!$OMP PARALLEL DO &
2446	!$OMP PRIVATE ( ij )
2447	DO ij = 1 , grid%num_tiles
2448	CALL set_physical_bc3d( grid%tke_1, 'p', config_flags, &
2449	ids, ide, jds, jde, kds, kde, &
2450	ims, ime, jms, jme, kms, kme, &
2451	ips, ipe, jps, jpe, kps, kpe, &
2452	grid%i_start(ij), grid%i_end(ij), &
2453	grid%j_start(ij), grid%j_end(ij), &
2454	k_start , k_end )
2455	END DO
2456	!$OMP END PARALLEL DO
2457
2458	!--- end of positive definite physics tendency update
2459
2460	END IF ! end if for tke_adv_opt
2461
2462	#ifdef DM_PARALLEL
2463	!
2464	! Stencils for patch communications (WCS, 29 June 2001)
2465	!
2466	! * * * * *
2467	! * * * * *
2468	! * * + * *
2469	! * * * * *
2470	! * * * * *
2471	!
2472	! ru_m x
2473	! rv_m x
2474	! ww_m x
2475	! mut x
2476	!
2477	!--------------------------------------------------------------
2478
2479	# include "HALO_EM_D.inc"
2480	! WCS addition 11/19/08
2481	# include "PERIOD_EM_DA.inc"
2482	#endif
2483
2484	!<DESCRIPTION>
2485	!<pre>
2486	! (4) Still within the RK loop, the scalar variables are advanced.
2487	!
2488	! For the moist and chem variables, each one is advanced
2489	! individually, using named loops "moist_variable_loop:"
2490	! and "chem_variable_loop:". Each RK substep begins by
2491	! calculating the advective tendency, and, for the first RK step,
2492	! 3D mixing (calling rk_scalar_tend) followed by an update
2493	! of the scalar (calling rk_update_scalar).
2494	!</pre>
2495	!</DESCRIPTION>
2496
2497	PRINT *,' Lluis PARAM_FIRST_SCALAR: ', PARAM_FIRST_SCALAR, &
2498	' p_qv: ',p_qv, ' p_qc: ',p_qc,' U(moist): ', UBOUND(moist)
2499	PRINT *,' moist: ',moist(config_flags%i_check_point, 2, &
2500	config_flags%j_check_point,:)
2501
2502	moist_scalar_advance: IF (num_3d_m >= PARAM_FIRST_SCALAR ) THEN
2503
2504	moist_variable_loop: DO im = PARAM_FIRST_SCALAR, num_3d_m
2505
2506	! adv_moist_cond is set in module_physics_init based on mp_physics choice
2507	! true except for Ferrier scheme
2508
2509	IF (grid%adv_moist_cond .or. im==p_qv ) THEN
2510
2511	!$OMP PARALLEL DO &
2512	!$OMP PRIVATE ( ij )
2513	moist_tile_loop_1: DO ij = 1 , grid%num_tiles
2514
2515	CALL wrf_debug ( 200 , ' call rk_scalar_tend' )
2516	tenddec = .false.
2517	PRINT *,' Lluis before rk_update_scalar moist im=',im, &
2518	' moist_old moist_tend moist________'
2519	DO iz = kms, kme
2520	PRINT *,moist_old(config_flags%i_check_point, iz, &
2521	config_flags%j_check_point,im), &
2522	moist_tend(config_flags%i_check_point, iz, &
2523	config_flags%j_check_point,im), &
2524	moist(config_flags%i_check_point, iz, &
2525	config_flags%j_check_point,im)
2526	END DO
2527
2528	BENCH_START(rk_scalar_tend_tim)
2529	CALL rk_scalar_tend ( im, im, config_flags, tenddec, &
2530	rk_step, dt_rk, &
2531	grid%ru_m, grid%rv_m, grid%ww_m, &
2532	grid%muts, grid%mub, grid%mu_1, &
2533	grid%alt, &
2534	moist_old(ims,kms,jms,im), &
2535	moist(ims,kms,jms,im), &
2536	moist_tend(ims,kms,jms,im), &
2537	advect_tend,h_tendency,z_tendency,grid%rqvften, &
2538	grid%qv_base, .true., grid%fnm, grid%fnp, &
2539	grid%msfux,grid%msfuy, grid%msfvx, grid%msfvx_inv,&
2540	grid%msfvy, grid%msftx,grid%msfty, &
2541	grid%rdx, grid%rdy, grid%rdn, grid%rdnw, grid%khdif, &
2542	grid%kvdif, grid%xkhh, &
2543	grid%diff_6th_opt, grid%diff_6th_factor, &
2544	config_flags%moist_adv_opt, &
2545	ids, ide, jds, jde, kds, kde, &
2546	ims, ime, jms, jme, kms, kme, &
2547	grid%i_start(ij), grid%i_end(ij), &
2548	grid%j_start(ij), grid%j_end(ij), &
2549	k_start , k_end )
2550
2551	BENCH_END(rk_scalar_tend_tim)
2552
2553	BENCH_START(rlx_bdy_scalar_tim)
2554	IF( ( config_flags%specified .or. config_flags%nested ) .and. rk_step == 1 ) THEN
2555	IF ( im .EQ. P_QV .OR. config_flags%nested ) THEN
2556	CALL relax_bdy_scalar ( moist_tend(ims,kms,jms,im), &
2557	moist(ims,kms,jms,im), grid%mut, &
2558	moist_bxs(jms,kms,1,im),moist_bxe(jms,kms,1,im), &
2559	moist_bys(ims,kms,1,im),moist_bye(ims,kms,1,im), &
2560	moist_btxs(jms,kms,1,im),moist_btxe(jms,kms,1,im), &
2561	moist_btys(ims,kms,1,im),moist_btye(ims,kms,1,im), &
2562	config_flags%spec_bdy_width, grid%spec_zone, grid%relax_zone, &
2563	grid%dtbc, grid%fcx, grid%gcx, &
2564	config_flags, &
2565	ids,ide, jds,jde, kds,kde, & ! domain dims
2566	ims,ime, jms,jme, kms,kme, & ! memory dims
2567	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
2568	grid%i_start(ij), grid%i_end(ij), &
2569	grid%j_start(ij), grid%j_end(ij), &
2570	k_start, k_end )
2571
2572	CALL spec_bdy_scalar ( moist_tend(ims,kms,jms,im), &
2573	moist_bxs(jms,kms,1,im),moist_bxe(jms,kms,1,im), &
2574	moist_bys(ims,kms,1,im),moist_bye(ims,kms,1,im), &
2575	moist_btxs(jms,kms,1,im),moist_btxe(jms,kms,1,im), &
2576	moist_btys(ims,kms,1,im),moist_btye(ims,kms,1,im), &
2577	config_flags%spec_bdy_width, grid%spec_zone, &
2578	config_flags, &
2579	ids,ide, jds,jde, kds,kde, & ! domain dims
2580	ims,ime, jms,jme, kms,kme, & ! memory dims
2581	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
2582	grid%i_start(ij), grid%i_end(ij), &
2583	grid%j_start(ij), grid%j_end(ij), &
2584	k_start, k_end )
2585	ENDIF
2586	ENDIF
2587	BENCH_END(rlx_bdy_scalar_tim)
2588
2589	ENDDO moist_tile_loop_1
2590	!$OMP END PARALLEL DO
2591
2592	!$OMP PARALLEL DO &
2593	!$OMP PRIVATE ( ij )
2594	moist_tile_loop_2: DO ij = 1 , grid%num_tiles
2595
2596	CALL wrf_debug ( 200 , ' call rk_update_scalar' )
2597	tenddec = .false.
2598
2599	BENCH_START(update_scal_tim)
2600	CALL rk_update_scalar( scs=im, sce=im, &
2601	scalar_1=moist_old(ims,kms,jms,im), &
2602	scalar_2=moist(ims,kms,jms,im), &
2603	sc_tend=moist_tend(ims,kms,jms,im), &
2604	! advh_t=advh_t(ims,kms,jms,1), &
2605	! advz_t=advz_t(ims,kms,jms,1), &
2606	advect_tend=advect_tend, &
2607	h_tendency=h_tendency, z_tendency=z_tendency, &
2608	msftx=grid%msftx,msfty=grid%msfty, &
2609	mu_old=grid%mu_1, mu_new=grid%mu_2, mu_base=grid%mub, &
2610	rk_step=rk_step, dt=dt_rk, spec_zone=grid%spec_zone, &
2611	config_flags=config_flags, tenddec=tenddec, &
2612	ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, &
2613	ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, &
2614	its=grid%i_start(ij), ite=grid%i_end(ij), &
2615	jts=grid%j_start(ij), jte=grid%j_end(ij), &
2616	kts=k_start , kte=k_end )
2617	BENCH_END(update_scal_tim)
2618	PRINT *,' Lluis after rk_update_scalar moist im=',im, &
2619	' moist_old moist_tend moist________'
2620	DO iz = kms, kme
2621	PRINT *,moist_old(config_flags%i_check_point, iz, &
2622	config_flags%j_check_point,im), &
2623	moist_tend(config_flags%i_check_point, iz, &
2624	config_flags%j_check_point,im), &
2625	moist(config_flags%i_check_point, iz, &
2626	config_flags%j_check_point,im)
2627	END DO
2628
2629	BENCH_START(flow_depbdy_tim)
2630	IF( config_flags%specified ) THEN
2631	IF(im .ne. P_QV)THEN
2632	CALL flow_dep_bdy ( moist(ims,kms,jms,im), &
2633	grid%ru_m, grid%rv_m, config_flags, &
2634	grid%spec_zone, &
2635	ids,ide, jds,jde, kds,kde, &
2636	ims,ime, jms,jme, kms,kme, &
2637	ips,ipe, jps,jpe, kps,kpe, &
2638	grid%i_start(ij), grid%i_end(ij), &
2639	grid%j_start(ij), grid%j_end(ij), &
2640	k_start, k_end )
2641	ENDIF
2642	ENDIF
2643	BENCH_END(flow_depbdy_tim)
2644
2645	ENDDO moist_tile_loop_2
2646	!$OMP END PARALLEL DO
2647
2648	ENDIF !-- if (grid%adv_moist_cond .or. im==p_qv ) then
2649
2650	ENDDO moist_variable_loop
2651
2652	ENDIF moist_scalar_advance
2653
2654	BENCH_START(tke_adv_tim)
2655	TKE_advance: IF (config_flags%km_opt .eq. 2) then
2656	#ifdef DM_PARALLEL
2657	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
2658	# include "HALO_EM_TKE_ADVECT_3.inc"
2659	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
2660	# include "HALO_EM_TKE_ADVECT_5.inc"
2661	ELSE
2662	WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order
2663	CALL wrf_error_fatal(TRIM(wrf_err_message))
2664	ENDIF
2665	#endif
2666	!$OMP PARALLEL DO &
2667	!$OMP PRIVATE ( ij )
2668	tke_tile_loop_1: DO ij = 1 , grid%num_tiles
2669
2670	CALL wrf_debug ( 200 , ' call rk_scalar_tend for tke' )
2671	tenddec = .false.
2672	CALL rk_scalar_tend ( 1, 1, config_flags, tenddec, &
2673	rk_step, dt_rk, &
2674	grid%ru_m, grid%rv_m, grid%ww_m, &
2675	grid%muts, grid%mub, grid%mu_1, &
2676	grid%alt, &
2677	grid%tke_1, &
2678	grid%tke_2, &
2679	tke_tend(ims,kms,jms), &
2680	advect_tend,h_tendency,z_tendency,grid%rqvften, &
2681	grid%qv_base, .false., grid%fnm, grid%fnp, &
2682	grid%msfux,grid%msfuy, grid%msfvx, grid%msfvx_inv, &
2683	grid%msfvy, grid%msftx,grid%msfty, &
2684	grid%rdx, grid%rdy, grid%rdn, grid%rdnw, grid%khdif, &
2685	grid%kvdif, grid%xkhh, &
2686	grid%diff_6th_opt, grid%diff_6th_factor, &
2687	config_flags%tke_adv_opt, &
2688	ids, ide, jds, jde, kds, kde, &
2689	ims, ime, jms, jme, kms, kme, &
2690	grid%i_start(ij), grid%i_end(ij), &
2691	grid%j_start(ij), grid%j_end(ij), &
2692	k_start , k_end )
2693
2694	ENDDO tke_tile_loop_1
2695	!$OMP END PARALLEL DO
2696
2697	!$OMP PARALLEL DO &
2698	!$OMP PRIVATE ( ij )
2699	tke_tile_loop_2: DO ij = 1 , grid%num_tiles
2700
2701	CALL wrf_debug ( 200 , ' call rk_update_scalar' )
2702	tenddec = .false.
2703	CALL rk_update_scalar( scs=1, sce=1, &
2704	scalar_1=grid%tke_1, &
2705	scalar_2=grid%tke_2, &
2706	sc_tend=tke_tend(ims,kms,jms), &
2707	! advh_t=advh_t(ims,kms,jms,1), &
2708	! advz_t=advz_t(ims,kms,jms,1), &
2709	advect_tend=advect_tend, &
2710	h_tendency=h_tendency, z_tendency=z_tendency, &
2711	msftx=grid%msftx,msfty=grid%msfty, &
2712	mu_old=grid%mu_1, mu_new=grid%mu_2, mu_base=grid%mub, &
2713	rk_step=rk_step, dt=dt_rk, spec_zone=grid%spec_zone, &
2714	config_flags=config_flags, tenddec=tenddec, &
2715	ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, &
2716	ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, &
2717	its=grid%i_start(ij), ite=grid%i_end(ij), &
2718	jts=grid%j_start(ij), jte=grid%j_end(ij), &
2719	kts=k_start , kte=k_end )
2720
2721	! bound the tke (greater than 0, less than tke_upper_bound)
2722
2723	CALL bound_tke( grid%tke_2, grid%tke_upper_bound, &
2724	ids, ide, jds, jde, kds, kde, &
2725	ims, ime, jms, jme, kms, kme, &
2726	grid%i_start(ij), grid%i_end(ij), &
2727	grid%j_start(ij), grid%j_end(ij), &
2728	k_start , k_end )
2729
2730	IF( config_flags%specified .or. config_flags%nested ) THEN
2731	CALL flow_dep_bdy ( grid%tke_2, &
2732	grid%ru_m, grid%rv_m, config_flags, &
2733	grid%spec_zone, &
2734	ids,ide, jds,jde, kds,kde, & ! domain dims
2735	ims,ime, jms,jme, kms,kme, & ! memory dims
2736	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
2737	grid%i_start(ij), grid%i_end(ij), &
2738	grid%j_start(ij), grid%j_end(ij), &
2739	k_start, k_end )
2740	ENDIF
2741	ENDDO tke_tile_loop_2
2742	!$OMP END PARALLEL DO
2743
2744	ENDIF TKE_advance
2745	BENCH_END(tke_adv_tim)
2746
2747	#ifdef WRF_CHEM
2748	! next the chemical species
2749	BENCH_START(chem_adv_tim)
2750	chem_scalar_advance: IF (num_3d_c >= PARAM_FIRST_SCALAR) THEN
2751
2752	chem_variable_loop: DO ic = PARAM_FIRST_SCALAR, num_3d_c
2753
2754	!$OMP PARALLEL DO &
2755	!$OMP PRIVATE ( ij )
2756	chem_tile_loop_1: DO ij = 1 , grid%num_tiles
2757
2758	CALL wrf_debug ( 200 , ' call rk_scalar_tend in chem_tile_loop_1' )
2759	tenddec = (( config_flags%chemdiag == USECHEMDIAG ) .and. &
2760	( adv_ct_indices(ic) >= PARAM_FIRST_SCALAR ))
2761	CALL rk_scalar_tend ( ic, ic, config_flags, tenddec, &
2762	rk_step, dt_rk, &
2763	grid%ru_m, grid%rv_m, grid%ww_m, &
2764	grid%muts, grid%mub, grid%mu_1, &
2765	grid%alt, &
2766	chem_old(ims,kms,jms,ic), &
2767	chem(ims,kms,jms,ic), &
2768	chem_tend(ims,kms,jms,ic), &
2769	advect_tend,h_tendency,z_tendency,grid%rqvften, &
2770	grid%qv_base, .false., grid%fnm, grid%fnp, &
2771	grid%msfux,grid%msfuy, grid%msfvx, grid%msfvx_inv, &
2772	grid%msfvy, grid%msftx,grid%msfty, &
2773	grid%rdx, grid%rdy, grid%rdn, grid%rdnw, &
2774	grid%khdif, grid%kvdif, grid%xkhh, &
2775	grid%diff_6th_opt, grid%diff_6th_factor, &
2776	config_flags%chem_adv_opt, &
2777	ids, ide, jds, jde, kds, kde, &
2778	ims, ime, jms, jme, kms, kme, &
2779	grid%i_start(ij), grid%i_end(ij), &
2780	grid%j_start(ij), grid%j_end(ij), &
2781	k_start , k_end )
2782	!
2783	! Currently, chemistry species with specified boundaries (i.e. the mother
2784	! domain) are being over written by flow_dep_bdy_chem. So, relax_bdy and
2785	! spec_bdy are only called for nests. For boundary conditions from global model or larger domain,
2786	! chem is uncoupled, and only used for one row/column on inflow (if have_bcs_chem=.true.)
2787	!
2788	IF( ( config_flags%nested ) .and. rk_step == 1 ) THEN
2789	IF(ic.eq.1)CALL wrf_debug ( 10 , ' have_bcs_chem' )
2790	CALL relax_bdy_scalar ( chem_tend(ims,kms,jms,ic), &
2791	chem(ims,kms,jms,ic), grid%mut, &
2792	chem_bxs(jms,kms,1,ic),chem_bxe(jms,kms,1,ic), &
2793	chem_bys(ims,kms,1,ic),chem_bye(ims,kms,1,ic), &
2794	chem_btxs(jms,kms,1,ic),chem_btxe(jms,kms,1,ic), &
2795	chem_btys(ims,kms,1,ic),chem_btye(ims,kms,1,ic), &
2796	config_flags%spec_bdy_width, grid%spec_zone, grid%relax_zone, &
2797	grid%dtbc, grid%fcx, grid%gcx, &
2798	config_flags, &
2799	ids,ide, jds,jde, kds,kde, &
2800	ims,ime, jms,jme, kms,kme, &
2801	ips,ipe, jps,jpe, kps,kpe, &
2802	grid%i_start(ij), grid%i_end(ij), &
2803	grid%j_start(ij), grid%j_end(ij), &
2804	k_start, k_end )
2805	CALL spec_bdy_scalar ( chem_tend(ims,kms,jms,ic), &
2806	chem_bxs(jms,kms,1,ic),chem_bxe(jms,kms,1,ic), &
2807	chem_bys(ims,kms,1,ic),chem_bye(ims,kms,1,ic), &
2808	chem_btxs(jms,kms,1,ic),chem_btxe(jms,kms,1,ic), &
2809	chem_btys(ims,kms,1,ic),chem_btye(ims,kms,1,ic), &
2810	config_flags%spec_bdy_width, grid%spec_zone, &
2811	config_flags, &
2812	ids,ide, jds,jde, kds,kde, &
2813	ims,ime, jms,jme, kms,kme, &
2814	ips,ipe, jps,jpe, kps,kpe, &
2815	grid%i_start(ij), grid%i_end(ij), &
2816	grid%j_start(ij), grid%j_end(ij), &
2817	k_start, k_end )
2818	ENDIF
2819
2820	ENDDO chem_tile_loop_1
2821	!$OMP END PARALLEL DO
2822
2823	!$OMP PARALLEL DO &
2824	!$OMP PRIVATE ( ij )
2825
2826	chem_tile_loop_2: DO ij = 1 , grid%num_tiles
2827
2828	CALL wrf_debug ( 200 , ' call rk_update_scalar' )
2829	tenddec = (( config_flags%chemdiag == USECHEMDIAG ) .and. &
2830	( adv_ct_indices(ic) >= PARAM_FIRST_SCALAR ))
2831	CALL rk_update_scalar( scs=ic, sce=ic, &
2832	scalar_1=chem_old(ims,kms,jms,ic), &
2833	scalar_2=chem(ims,kms,jms,ic), &
2834	sc_tend=chem_tend(ims,kms,jms,ic), &
2835	advh_t=advh_ct(ims,kms,jms,adv_ct_indices(ic)), &
2836	advz_t=advz_ct(ims,kms,jms,adv_ct_indices(ic)), &
2837	advect_tend=advect_tend, &
2838	h_tendency=h_tendency, z_tendency=z_tendency, &
2839	msftx=grid%msftx,msfty=grid%msfty, &
2840	mu_old=grid%mu_1, mu_new=grid%mu_2, mu_base=grid%mub, &
2841	rk_step=rk_step, dt=dt_rk, spec_zone=grid%spec_zone, &
2842	config_flags=config_flags, tenddec=tenddec, &
2843	ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, &
2844	ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, &
2845	its=grid%i_start(ij), ite=grid%i_end(ij), &
2846	jts=grid%j_start(ij), jte=grid%j_end(ij), &
2847	kts=k_start , kte=k_end )
2848
2849	IF( config_flags%specified ) THEN
2850	CALL flow_dep_bdy_chem( chem(ims,kms,jms,ic), &
2851	chem_bxs(jms,kms,1,ic), chem_btxs(jms,kms,1,ic), &
2852	chem_bxe(jms,kms,1,ic), chem_btxe(jms,kms,1,ic), &
2853	chem_bys(ims,kms,1,ic), chem_btys(ims,kms,1,ic), &
2854	chem_bye(ims,kms,1,ic), chem_btye(ims,kms,1,ic), &
2855	dt_rk+grid%dtbc, &
2856	config_flags%spec_bdy_width,grid%z, &
2857	grid%have_bcs_chem, &
2858	grid%ru_m, grid%rv_m, config_flags,grid%alt, &
2859	grid%t_1,grid%pb,grid%p,t0,p1000mb,rcp,grid%ph_2,grid%phb,g, &
2860	grid%spec_zone,ic, &
2861	ids,ide, jds,jde, kds,kde, & ! domain dims
2862	ims,ime, jms,jme, kms,kme, & ! memory dims
2863	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
2864	grid%i_start(ij), grid%i_end(ij), &
2865	grid%j_start(ij), grid%j_end(ij), &
2866	k_start, k_end )
2867	ENDIF
2868	ENDDO chem_tile_loop_2
2869	!$OMP END PARALLEL DO
2870
2871	ENDDO chem_variable_loop
2872	ENDIF chem_scalar_advance
2873	BENCH_END(chem_adv_tim)
2874	#endif
2875	! next the chemical species
2876	BENCH_START(tracer_adv_tim)
2877	tracer_advance: IF (num_tracer >= PARAM_FIRST_SCALAR) THEN
2878
2879	tracer_variable_loop: DO ic = PARAM_FIRST_SCALAR, num_tracer
2880
2881	!$OMP PARALLEL DO &
2882	!$OMP PRIVATE ( ij )
2883	tracer_tile_loop_1: DO ij = 1 , grid%num_tiles
2884
2885	CALL wrf_debug ( 15 , ' call rk_scalar_tend in tracer_tile_loop_1' )
2886	tenddec = .false.
2887	CALL rk_scalar_tend ( ic, ic, config_flags, tenddec, &
2888	rk_step, dt_rk, &
2889	grid%ru_m, grid%rv_m, grid%ww_m, &
2890	grid%muts, grid%mub, grid%mu_1, &
2891	grid%alt, &
2892	tracer_old(ims,kms,jms,ic), &
2893	tracer(ims,kms,jms,ic), &
2894	tracer_tend(ims,kms,jms,ic), &
2895	advect_tend,h_tendency,z_tendency,grid%rqvften, &
2896	grid%qv_base, .false., grid%fnm, grid%fnp, &
2897	grid%msfux,grid%msfuy, grid%msfvx, grid%msfvx_inv, &
2898	grid%msfvy, grid%msftx,grid%msfty, &
2899	grid%rdx, grid%rdy, grid%rdn, grid%rdnw, &
2900	grid%khdif, grid%kvdif, grid%xkhh, &
2901	grid%diff_6th_opt, grid%diff_6th_factor, &
2902	config_flags%tracer_adv_opt, &
2903	ids, ide, jds, jde, kds, kde, &
2904	ims, ime, jms, jme, kms, kme, &
2905	grid%i_start(ij), grid%i_end(ij), &
2906	grid%j_start(ij), grid%j_end(ij), &
2907	k_start , k_end )
2908	!
2909	! Currently, chemistry species with specified boundaries (i.e. the mother
2910	! domain) are being over written by flow_dep_bdy_chem. So, relax_bdy and
2911	! spec_bdy are only called for nests. For boundary conditions from global model or larger domain,
2912	! chem is uncoupled, and only used for one row/column on inflow (if have_bcs_chem=.true.)
2913	!
2914	IF( ( config_flags%nested ) .and. rk_step == 1 ) THEN
2915	IF(ic.eq.1)CALL wrf_debug ( 10 , ' have_bcs_tracer' )
2916	CALL relax_bdy_scalar ( tracer_tend(ims,kms,jms,ic), &
2917	tracer(ims,kms,jms,ic), grid%mut, &
2918	tracer_bxs(jms,kms,1,ic),tracer_bxe(jms,kms,1,ic), &
2919	tracer_bys(ims,kms,1,ic),tracer_bye(ims,kms,1,ic), &
2920	tracer_btxs(jms,kms,1,ic),tracer_btxe(jms,kms,1,ic), &
2921	tracer_btys(ims,kms,1,ic),tracer_btye(ims,kms,1,ic), &
2922	config_flags%spec_bdy_width, grid%spec_zone, grid%relax_zone, &
2923	grid%dtbc, grid%fcx, grid%gcx, &
2924	config_flags, &
2925	ids,ide, jds,jde, kds,kde, &
2926	ims,ime, jms,jme, kms,kme, &
2927	ips,ipe, jps,jpe, kps,kpe, &
2928	grid%i_start(ij), grid%i_end(ij), &
2929	grid%j_start(ij), grid%j_end(ij), &
2930	k_start, k_end )
2931	CALL spec_bdy_scalar ( tracer_tend(ims,kms,jms,ic), &
2932	tracer_bxs(jms,kms,1,ic),tracer_bxe(jms,kms,1,ic), &
2933	tracer_bys(ims,kms,1,ic),tracer_bye(ims,kms,1,ic), &
2934	tracer_btxs(jms,kms,1,ic),tracer_btxe(jms,kms,1,ic), &
2935	tracer_btys(ims,kms,1,ic),tracer_btye(ims,kms,1,ic), &
2936	config_flags%spec_bdy_width, grid%spec_zone, &
2937	config_flags, &
2938	ids,ide, jds,jde, kds,kde, &
2939	ims,ime, jms,jme, kms,kme, &
2940	ips,ipe, jps,jpe, kps,kpe, &
2941	grid%i_start(ij), grid%i_end(ij), &
2942	grid%j_start(ij), grid%j_end(ij), &
2943	k_start, k_end )
2944	ENDIF
2945
2946	ENDDO tracer_tile_loop_1
2947	!$OMP END PARALLEL DO
2948
2949	!$OMP PARALLEL DO &
2950	!$OMP PRIVATE ( ij )
2951
2952	tracer_tile_loop_2: DO ij = 1 , grid%num_tiles
2953
2954	CALL wrf_debug ( 200 , ' call rk_update_scalar' )
2955	tenddec = .false.
2956	CALL rk_update_scalar( scs=ic, sce=ic, &
2957	scalar_1=tracer_old(ims,kms,jms,ic), &
2958	scalar_2=tracer(ims,kms,jms,ic), &
2959	sc_tend=tracer_tend(ims,kms,jms,ic), &
2960	! advh_t=advh_t(ims,kms,jms,1), &
2961	! advz_t=advz_t(ims,kms,jms,1), &
2962	advect_tend=advect_tend, &
2963	h_tendency=h_tendency, z_tendency=z_tendency, &
2964	msftx=grid%msftx,msfty=grid%msfty, &
2965	mu_old=grid%mu_1, mu_new=grid%mu_2, mu_base=grid%mub, &
2966	rk_step=rk_step, dt=dt_rk, spec_zone=grid%spec_zone, &
2967	config_flags=config_flags, tenddec=tenddec, &
2968	ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, &
2969	ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, &
2970	its=grid%i_start(ij), ite=grid%i_end(ij), &
2971	jts=grid%j_start(ij), jte=grid%j_end(ij), &
2972	kts=k_start , kte=k_end )
2973
2974	IF( config_flags%specified ) THEN
2975	#ifdef WRF_CHEM
2976	CALL flow_dep_bdy_tracer( tracer(ims,kms,jms,ic), &
2977	tracer_bxs(jms,kms,1,ic), tracer_btxs(jms,kms,1,ic), &
2978	tracer_bxe(jms,kms,1,ic), tracer_btxe(jms,kms,1,ic), &
2979	tracer_bys(ims,kms,1,ic), tracer_btys(ims,kms,1,ic), &
2980	tracer_bye(ims,kms,1,ic), tracer_btye(ims,kms,1,ic), &
2981	dt_rk+grid%dtbc, &
2982	config_flags%spec_bdy_width,grid%z, &
2983	grid%have_bcs_tracer, &
2984	grid%ru_m, grid%rv_m, config_flags%tracer_opt,grid%alt, &
2985	grid%t_1,grid%pb,grid%p,t0,p1000mb,rcp,grid%ph_2,grid%phb,g, &
2986	grid%spec_zone,ic, &
2987	ids,ide, jds,jde, kds,kde, & ! domain dims
2988	ims,ime, jms,jme, kms,kme, & ! memory dims
2989	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
2990	grid%i_start(ij), grid%i_end(ij), &
2991	grid%j_start(ij), grid%j_end(ij), &
2992	k_start, k_end )
2993	#else
2994	CALL flow_dep_bdy ( tracer(ims,kms,jms,ic), &
2995	grid%ru_m, grid%rv_m, config_flags, &
2996	grid%spec_zone, &
2997	ids,ide, jds,jde, kds,kde, & ! domain dims
2998	ims,ime, jms,jme, kms,kme, & ! memory dims
2999	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
3000	grid%i_start(ij), grid%i_end(ij), &
3001	grid%j_start(ij), grid%j_end(ij), &
3002	k_start, k_end )
3003	#endif
3004	ENDIF
3005	ENDDO tracer_tile_loop_2
3006	!$OMP END PARALLEL DO
3007
3008	ENDDO tracer_variable_loop
3009	ENDIF tracer_advance
3010	BENCH_END(tracer_adv_tim)
3011
3012	! next the other scalar species
3013	other_scalar_advance: IF (num_3d_s >= PARAM_FIRST_SCALAR) THEN
3014
3015	scalar_variable_loop: do is = PARAM_FIRST_SCALAR, num_3d_s
3016	!$OMP PARALLEL DO &
3017	!$OMP PRIVATE ( ij )
3018	scalar_tile_loop_1: DO ij = 1 , grid%num_tiles
3019
3020	CALL wrf_debug ( 200 , ' call rk_scalar_tend' )
3021	tenddec = .false.
3022	CALL rk_scalar_tend ( is, is, config_flags, tenddec, &
3023	rk_step, dt_rk, &
3024	grid%ru_m, grid%rv_m, grid%ww_m, &
3025	grid%muts, grid%mub, grid%mu_1, &
3026	grid%alt, &
3027	scalar_old(ims,kms,jms,is), &
3028	scalar(ims,kms,jms,is), &
3029	scalar_tend(ims,kms,jms,is), &
3030	advect_tend,h_tendency,z_tendency,grid%rqvften, &
3031	grid%qv_base, .false., grid%fnm, grid%fnp, &
3032	grid%msfux,grid%msfuy, grid%msfvx, grid%msfvx_inv, &
3033	grid%msfvy, grid%msftx,grid%msfty, &
3034	grid%rdx, grid%rdy, grid%rdn, grid%rdnw, &
3035	grid%khdif, grid%kvdif, grid%xkhh, &
3036	grid%diff_6th_opt, grid%diff_6th_factor, &
3037	config_flags%scalar_adv_opt, &
3038	ids, ide, jds, jde, kds, kde, &
3039	ims, ime, jms, jme, kms, kme, &
3040	grid%i_start(ij), grid%i_end(ij), &
3041	grid%j_start(ij), grid%j_end(ij), &
3042	k_start , k_end )
3043
3044	IF( config_flags%nested .and. (rk_step == 1) ) THEN
3045
3046	CALL relax_bdy_scalar ( scalar_tend(ims,kms,jms,is), &
3047	scalar(ims,kms,jms,is), grid%mut, &
3048	scalar_bxs(jms,kms,1,is),scalar_bxe(jms,kms,1,is), &
3049	scalar_bys(ims,kms,1,is),scalar_bye(ims,kms,1,is), &
3050	scalar_btxs(jms,kms,1,is),scalar_btxe(jms,kms,1,is), &
3051	scalar_btys(ims,kms,1,is),scalar_btye(ims,kms,1,is), &
3052	config_flags%spec_bdy_width, grid%spec_zone, grid%relax_zone, &
3053	grid%dtbc, grid%fcx, grid%gcx, &
3054	config_flags, &
3055	ids,ide, jds,jde, kds,kde, &
3056	ims,ime, jms,jme, kms,kme, &
3057	ips,ipe, jps,jpe, kps,kpe, &
3058	grid%i_start(ij), grid%i_end(ij), &
3059	grid%j_start(ij), grid%j_end(ij), &
3060	k_start, k_end )
3061
3062	CALL spec_bdy_scalar ( scalar_tend(ims,kms,jms,is), &
3063	scalar_bxs(jms,kms,1,is),scalar_bxe(jms,kms,1,is), &
3064	scalar_bys(ims,kms,1,is),scalar_bye(ims,kms,1,is), &
3065	scalar_btxs(jms,kms,1,is),scalar_btxe(jms,kms,1,is), &
3066	scalar_btys(ims,kms,1,is),scalar_btye(ims,kms,1,is), &
3067	config_flags%spec_bdy_width, grid%spec_zone, &
3068	config_flags, &
3069	ids,ide, jds,jde, kds,kde, &
3070	ims,ime, jms,jme, kms,kme, &
3071	ips,ipe, jps,jpe, kps,kpe, &
3072	grid%i_start(ij), grid%i_end(ij), &
3073	grid%j_start(ij), grid%j_end(ij), &
3074	k_start, k_end )
3075
3076	ENDIF ! b.c test for chem nested boundary condition
3077
3078	ENDDO scalar_tile_loop_1
3079	!$OMP END PARALLEL DO
3080
3081	!$OMP PARALLEL DO &
3082	!$OMP PRIVATE ( ij )
3083	scalar_tile_loop_2: DO ij = 1 , grid%num_tiles
3084
3085	CALL wrf_debug ( 200 , ' call rk_update_scalar' )
3086	tenddec = .false.
3087	CALL rk_update_scalar( scs=is, sce=is, &
3088	scalar_1=scalar_old(ims,kms,jms,is), &
3089	scalar_2=scalar(ims,kms,jms,is), &
3090	sc_tend=scalar_tend(ims,kms,jms,is), &
3091	! advh_t=advh_t(ims,kms,jms,1), &
3092	! advz_t=advz_t(ims,kms,jms,1), &
3093	advect_tend=advect_tend, &
3094	h_tendency=h_tendency, z_tendency=z_tendency, &
3095	msftx=grid%msftx,msfty=grid%msfty, &
3096	mu_old=grid%mu_1, mu_new=grid%mu_2, mu_base=grid%mub, &
3097	rk_step=rk_step, dt=dt_rk, spec_zone=grid%spec_zone, &
3098	config_flags=config_flags, tenddec=tenddec, &
3099	ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, &
3100	ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, &
3101	its=grid%i_start(ij), ite=grid%i_end(ij), &
3102	jts=grid%j_start(ij), jte=grid%j_end(ij), &
3103	kts=k_start , kte=k_end )
3104
3105	IF( config_flags%specified ) THEN
3106
3107	CALL flow_dep_bdy ( scalar(ims,kms,jms,is), &
3108	grid%ru_m, grid%rv_m, config_flags, &
3109	grid%spec_zone, &
3110	ids,ide, jds,jde, kds,kde, & ! domain dims
3111	ims,ime, jms,jme, kms,kme, & ! memory dims
3112	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
3113	grid%i_start(ij), grid%i_end(ij), &
3114	grid%j_start(ij), grid%j_end(ij), &
3115	k_start, k_end )
3116	ENDIF
3117
3118	ENDDO scalar_tile_loop_2
3119	!$OMP END PARALLEL DO
3120
3121	ENDDO scalar_variable_loop
3122
3123	ENDIF other_scalar_advance
3124
3125	! update the pressure and density at the new time level
3126
3127	!$OMP PARALLEL DO &
3128	!$OMP PRIVATE ( ij )
3129	DO ij = 1 , grid%num_tiles
3130
3131	BENCH_START(calc_p_rho_tim)
3132
3133	CALL calc_p_rho_phi( moist, num_3d_m, &
3134	grid%al, grid%alb, grid%mu_2, grid%muts, &
3135	grid%ph_2, grid%p, grid%pb, grid%t_2, &
3136	p0, t0, grid%znu, grid%dnw, grid%rdnw, &
3137	grid%rdn, config_flags%non_hydrostatic, &
3138	ids, ide, jds, jde, kds, kde, &
3139	ims, ime, jms, jme, kms, kme, &
3140	grid%i_start(ij), grid%i_end(ij), &
3141	grid%j_start(ij), grid%j_end(ij), &
3142	k_start , k_end )
3143
3144	BENCH_END(calc_p_rho_tim)
3145
3146	ENDDO
3147	!$OMP END PARALLEL DO
3148
3149	! Reset the boundary conditions if there is another corrector step.
3150	! (rk_step < rk_order), else we'll handle it at the end of everything
3151	! (after the split physics, before exiting the timestep).
3152
3153	rk_step_1_check: IF ( rk_step < rk_order ) THEN
3154
3155	!-----------------------------------------------------------
3156	! rk3 substep polar filter for scalars (moist,chem,scalar)
3157	!-----------------------------------------------------------
3158
3159	IF (config_flags%polar) THEN
3160	IF ( num_3d_m >= PARAM_FIRST_SCALAR ) THEN
3161	CALL wrf_debug ( 200 , ' call filter moist ' )
3162	DO im = PARAM_FIRST_SCALAR, num_3d_m
3163	CALL couple_scalars_for_filter ( FIELD=moist(ims,kms,jms,im) &
3164	,MU=grid%mu_2 , MUB=grid%mub &
3165	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3166	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3167	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe )
3168	CALL pxft ( grid=grid &
3169	,lineno=__LINE__ &
3170	,flag_uv = 0 &
3171	,flag_rurv = 0 &
3172	,flag_wph = 0 &
3173	,flag_ww = 0 &
3174	,flag_t = 0 &
3175	,flag_mu = 0 &
3176	,flag_mut = 0 &
3177	,flag_moist = im &
3178	,flag_chem = 0 &
3179	,flag_scalar = 0 &
3180	,flag_tracer = 0 &
3181	,positive_definite=.FALSE. &
3182	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
3183	,fft_filter_lat = config_flags%fft_filter_lat &
3184	,dclat = dclat &
3185	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3186	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3187	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
3188	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
3189	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
3190	CALL uncouple_scalars_for_filter ( FIELD=moist(ims,kms,jms,im) &
3191	,MU=grid%mu_2 , MUB=grid%mub &
3192	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3193	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3194	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe )
3195	END DO
3196	END IF
3197
3198	IF ( num_3d_c >= PARAM_FIRST_SCALAR ) THEN
3199	CALL wrf_debug ( 200 , ' call filter chem ' )
3200	DO im = PARAM_FIRST_SCALAR, num_3d_c
3201	CALL couple_scalars_for_filter ( FIELD=chem(ims,kms,jms,im) &
3202	,MU=grid%mu_2 , MUB=grid%mub &
3203	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3204	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3205	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe )
3206	CALL pxft ( grid=grid &
3207	,lineno=__LINE__ &
3208	,flag_uv = 0 &
3209	,flag_rurv = 0 &
3210	,flag_wph = 0 &
3211	,flag_ww = 0 &
3212	,flag_t = 0 &
3213	,flag_mu = 0 &
3214	,flag_mut = 0 &
3215	,flag_moist = 0 &
3216	,flag_chem = im &
3217	,flag_tracer = 0 &
3218	,flag_scalar = 0 &
3219	,positive_definite=.FALSE. &
3220	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
3221	,fft_filter_lat = config_flags%fft_filter_lat &
3222	,dclat = dclat &
3223	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3224	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3225	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
3226	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
3227	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
3228	CALL uncouple_scalars_for_filter ( FIELD=chem(ims,kms,jms,im) &
3229	,MU=grid%mu_2 , MUB=grid%mub &
3230	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3231	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3232	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe )
3233	END DO
3234	END IF
3235	IF ( num_tracer >= PARAM_FIRST_SCALAR ) THEN
3236	CALL wrf_debug ( 200 , ' call filter tracer ' )
3237	DO im = PARAM_FIRST_SCALAR, num_tracer
3238	CALL couple_scalars_for_filter ( FIELD=tracer(ims,kms,jms,im) &
3239	,MU=grid%mu_2 , MUB=grid%mub &
3240	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3241	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3242	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe )
3243	CALL pxft ( grid=grid &
3244	,lineno=__LINE__ &
3245	,flag_uv = 0 &
3246	,flag_rurv = 0 &
3247	,flag_wph = 0 &
3248	,flag_ww = 0 &
3249	,flag_t = 0 &
3250	,flag_mu = 0 &
3251	,flag_mut = 0 &
3252	,flag_moist = 0 &
3253	,flag_chem = 0 &
3254	,flag_tracer = im &
3255	,flag_scalar = 0 &
3256	,positive_definite=.FALSE. &
3257	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
3258	,fft_filter_lat = config_flags%fft_filter_lat &
3259	,dclat = dclat &
3260	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3261	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3262	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
3263	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
3264	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
3265	CALL uncouple_scalars_for_filter ( FIELD=tracer(ims,kms,jms,im) &
3266	,MU=grid%mu_2 , MUB=grid%mub &
3267	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3268	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3269	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe )
3270	END DO
3271	END IF
3272
3273	IF ( num_3d_s >= PARAM_FIRST_SCALAR ) THEN
3274	CALL wrf_debug ( 200 , ' call filter scalar ' )
3275	DO im = PARAM_FIRST_SCALAR, num_3d_s
3276	CALL couple_scalars_for_filter ( FIELD=scalar(ims,kms,jms,im) &
3277	,MU=grid%mu_2 , MUB=grid%mub &
3278	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3279	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3280	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe )
3281	CALL pxft ( grid=grid &
3282	,lineno=__LINE__ &
3283	,flag_uv = 0 &
3284	,flag_rurv = 0 &
3285	,flag_wph = 0 &
3286	,flag_ww = 0 &
3287	,flag_t = 0 &
3288	,flag_mu = 0 &
3289	,flag_mut = 0 &
3290	,flag_moist = 0 &
3291	,flag_chem = 0 &
3292	,flag_tracer = 0 &
3293	,flag_scalar = im &
3294	,positive_definite=.FALSE. &
3295	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
3296	,fft_filter_lat = config_flags%fft_filter_lat &
3297	,dclat = dclat &
3298	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3299	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3300	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
3301	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
3302	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
3303	CALL uncouple_scalars_for_filter ( FIELD=scalar(ims,kms,jms,im) &
3304	,MU=grid%mu_2 , MUB=grid%mub &
3305	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3306	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3307	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe )
3308	END DO
3309	END IF
3310	END IF ! polar filter test
3311
3312	!-----------------------------------------------------------
3313	! END rk3 substep polar filter for scalars (moist,chem,scalar)
3314	!-----------------------------------------------------------
3315
3316	!-----------------------------------------------------------
3317	! Stencils for patch communications (WCS, 29 June 2001)
3318	!
3319	! here's where we need a wide comm stencil - these are the
3320	! uncoupled variables so are used for high order calc in
3321	! advection and mixong routines.
3322	!
3323	!
3324	! * * * * * * *
3325	! * * * * * * * * * * * *
3326	! * * * * * * * * * * * * *
3327	! * + * * * + * * * * * + * * *
3328	! * * * * * * * * * * * * *
3329	! * * * * * * * * * * * *
3330	! * * * * * * *
3331	!
3332	! al x
3333	!
3334	! 2D variable
3335	! mu_2 x
3336	!
3337	! (adv order <=4)
3338	! u_2 x
3339	! v_2 x
3340	! w_2 x
3341	! t_2 x
3342	! ph_2 x
3343	!
3344	! (adv order <=6)
3345	! u_2 x
3346	! v_2 x
3347	! w_2 x
3348	! t_2 x
3349	! ph_2 x
3350	!
3351	! 4D variable
3352	! moist x
3353	! chem x
3354	! scalar x
3355
3356	#ifdef DM_PARALLEL
3357	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
3358	# include "HALO_EM_D2_3.inc"
3359	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
3360	# include "HALO_EM_D2_5.inc"
3361	ELSE
3362	WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order
3363	CALL wrf_error_fatal(TRIM(wrf_err_message))
3364	ENDIF
3365	# include "PERIOD_BDY_EM_D.inc"
3366	# include "PERIOD_BDY_EM_MOIST2.inc"
3367	# include "PERIOD_BDY_EM_CHEM2.inc"
3368	# include "PERIOD_BDY_EM_TRACER2.inc"
3369	# include "PERIOD_BDY_EM_SCALAR2.inc"
3370	#endif
3371
3372	BENCH_START(bc_end_tim)
3373	!$OMP PARALLEL DO &
3374	!$OMP PRIVATE ( ij )
3375	tile_bc_loop_1: DO ij = 1 , grid%num_tiles
3376	CALL wrf_debug ( 200 , ' call rk_phys_bc_dry_2' )
3377
3378	CALL rk_phys_bc_dry_2( config_flags, &
3379	grid%u_2, grid%v_2, grid%w_2, &
3380	grid%t_2, grid%ph_2, grid%mu_2, &
3381	ids, ide, jds, jde, kds, kde, &
3382	ims, ime, jms, jme, kms, kme, &
3383	ips, ipe, jps, jpe, kps, kpe, &
3384	grid%i_start(ij), grid%i_end(ij), &
3385	grid%j_start(ij), grid%j_end(ij), &
3386	k_start , k_end )
3387
3388	BENCH_START(diag_w_tim)
3389	IF (.not. config_flags%non_hydrostatic) THEN
3390	CALL diagnose_w( ph_tend, grid%ph_2, grid%ph_1, grid%w_2, grid%muts, dt_rk, &
3391	grid%u_2, grid%v_2, grid%ht, &
3392	grid%cf1, grid%cf2, grid%cf3, grid%rdx, grid%rdy, grid%msftx, grid%msfty, &
3393	ids, ide, jds, jde, kds, kde, &
3394	ims, ime, jms, jme, kms, kme, &
3395	grid%i_start(ij), grid%i_end(ij), &
3396	grid%j_start(ij), grid%j_end(ij), &
3397	k_start , k_end )
3398	ENDIF
3399	BENCH_END(diag_w_tim)
3400
3401	IF (num_3d_m >= PARAM_FIRST_SCALAR) THEN
3402
3403	moisture_loop_bdy_1 : DO im = PARAM_FIRST_SCALAR , num_3d_m
3404
3405	CALL set_physical_bc3d( moist(ims,kms,jms,im), 'p', config_flags, &
3406	ids, ide, jds, jde, kds, kde, &
3407	ims, ime, jms, jme, kms, kme, &
3408	ips, ipe, jps, jpe, kps, kpe, &
3409	grid%i_start(ij), grid%i_end(ij), &
3410	grid%j_start(ij), grid%j_end(ij), &
3411	k_start , k_end )
3412	END DO moisture_loop_bdy_1
3413
3414	ENDIF
3415
3416	IF (num_3d_c >= PARAM_FIRST_SCALAR) THEN
3417
3418	chem_species_bdy_loop_1 : DO ic = PARAM_FIRST_SCALAR , num_3d_c
3419
3420	CALL set_physical_bc3d( chem(ims,kms,jms,ic), 'p', config_flags, &
3421	ids, ide, jds, jde, kds, kde, &
3422	ims, ime, jms, jme, kms, kme, &
3423	ips, ipe, jps, jpe, kps, kpe, &
3424	grid%i_start(ij), grid%i_end(ij), &
3425	grid%j_start(ij), grid%j_end(ij), &
3426	k_start , k_end-1 )
3427
3428	END DO chem_species_bdy_loop_1
3429
3430	END IF
3431
3432	IF (num_tracer >= PARAM_FIRST_SCALAR) THEN
3433
3434	tracer_species_bdy_loop_1 : DO ic = PARAM_FIRST_SCALAR , num_tracer
3435
3436	CALL set_physical_bc3d( tracer(ims,kms,jms,ic), 'p', config_flags, &
3437	ids, ide, jds, jde, kds, kde, &
3438	ims, ime, jms, jme, kms, kme, &
3439	ips, ipe, jps, jpe, kps, kpe, &
3440	grid%i_start(ij), grid%i_end(ij), &
3441	grid%j_start(ij), grid%j_end(ij), &
3442	k_start , k_end-1 )
3443
3444	END DO tracer_species_bdy_loop_1
3445
3446	END IF
3447
3448	IF (num_3d_s >= PARAM_FIRST_SCALAR) THEN
3449
3450	scalar_species_bdy_loop_1 : DO is = PARAM_FIRST_SCALAR , num_3d_s
3451
3452	CALL set_physical_bc3d( scalar(ims,kms,jms,is), 'p', config_flags, &
3453	ids, ide, jds, jde, kds, kde, &
3454	ims, ime, jms, jme, kms, kme, &
3455	ips, ipe, jps, jpe, kps, kpe, &
3456	grid%i_start(ij), grid%i_end(ij), &
3457	grid%j_start(ij), grid%j_end(ij), &
3458	k_start , k_end-1 )
3459
3460	END DO scalar_species_bdy_loop_1
3461
3462	END IF
3463
3464	IF (config_flags%km_opt .eq. 2) THEN
3465
3466	CALL set_physical_bc3d( grid%tke_2 , 'p', config_flags, &
3467	ids, ide, jds, jde, kds, kde, &
3468	ims, ime, jms, jme, kms, kme, &
3469	ips, ipe, jps, jpe, kps, kpe, &
3470	grid%i_start(ij), grid%i_end(ij), &
3471	grid%j_start(ij), grid%j_end(ij), &
3472	k_start , k_end )
3473	END IF
3474
3475	END DO tile_bc_loop_1
3476	!$OMP END PARALLEL DO
3477	BENCH_END(bc_end_tim)
3478
3479
3480	#ifdef DM_PARALLEL
3481
3482	! * * * * *
3483	! * * * * * * * * *
3484	! * + * * + * * * + * *
3485	! * * * * * * * * *
3486	! * * * * *
3487
3488	! moist, chem, scalar, tke x
3489
3490
3491	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
3492	IF ( (config_flags%tke_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3493	# include "HALO_EM_TKE_5.inc"
3494	ELSE
3495	# include "HALO_EM_TKE_3.inc"
3496	ENDIF
3497	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
3498	IF ( (config_flags%tke_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3499	# include "HALO_EM_TKE_7.inc"
3500	ELSE
3501	# include "HALO_EM_TKE_5.inc"
3502	ENDIF
3503	ELSE
3504	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
3505	CALL wrf_error_fatal(TRIM(wrf_err_message))
3506	ENDIF
3507
3508	IF ( num_moist .GE. PARAM_FIRST_SCALAR ) THEN
3509	IF ( config_flags%h_sca_adv_order <= 4 ) THEN
3510	IF ( (config_flags%moist_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3511	# include "HALO_EM_MOIST_E_5.inc"
3512	ELSE
3513	# include "HALO_EM_MOIST_E_3.inc"
3514	END IF
3515	ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN
3516	IF ( (config_flags%moist_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3517	# include "HALO_EM_MOIST_E_7.inc"
3518	ELSE
3519	# include "HALO_EM_MOIST_E_5.inc"
3520	END IF
3521	ELSE
3522	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
3523	CALL wrf_error_fatal(TRIM(wrf_err_message))
3524	ENDIF
3525	ENDIF
3526	IF ( num_chem >= PARAM_FIRST_SCALAR ) THEN
3527	IF ( config_flags%h_sca_adv_order <= 4 ) THEN
3528	IF ( (config_flags%chem_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3529	# include "HALO_EM_CHEM_E_5.inc"
3530	ELSE
3531	# include "HALO_EM_CHEM_E_3.inc"
3532	ENDIF
3533	ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN
3534	IF ( (config_flags%chem_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3535	# include "HALO_EM_CHEM_E_7.inc"
3536	ELSE
3537	# include "HALO_EM_CHEM_E_5.inc"
3538	ENDIF
3539	ELSE
3540	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
3541	CALL wrf_error_fatal(TRIM(wrf_err_message))
3542	ENDIF
3543	ENDIF
3544	IF ( num_tracer >= PARAM_FIRST_SCALAR ) THEN
3545	IF ( config_flags%h_sca_adv_order <= 4 ) THEN
3546	IF ( (config_flags%tracer_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3547	# include "HALO_EM_TRACER_E_5.inc"
3548	ELSE
3549	# include "HALO_EM_TRACER_E_3.inc"
3550	ENDIF
3551	ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN
3552	IF ( (config_flags%tracer_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3553	# include "HALO_EM_TRACER_E_7.inc"
3554	ELSE
3555	# include "HALO_EM_TRACER_E_5.inc"
3556	ENDIF
3557	ELSE
3558	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
3559	CALL wrf_error_fatal(TRIM(wrf_err_message))
3560	ENDIF
3561	ENDIF
3562	IF ( num_scalar >= PARAM_FIRST_SCALAR ) THEN
3563	IF ( config_flags%h_sca_adv_order <= 4 ) THEN
3564	IF ( (config_flags%scalar_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3565	# include "HALO_EM_SCALAR_E_5.inc"
3566	ELSE
3567	# include "HALO_EM_SCALAR_E_3.inc"
3568	ENDIF
3569	ELSE IF ( config_flags%h_sca_adv_order <= 6 ) THEN
3570	IF ( (config_flags%scalar_adv_opt /= ORIGINAL) .and. (rk_step == rk_order-1) ) THEN
3571	# include "HALO_EM_SCALAR_E_7.inc"
3572	ELSE
3573	# include "HALO_EM_SCALAR_E_5.inc"
3574	ENDIF
3575	ELSE
3576	WRITE(wrf_err_message,*)'solve_em: invalid h_sca_adv_order = ',config_flags%h_sca_adv_order
3577	CALL wrf_error_fatal(TRIM(wrf_err_message))
3578	ENDIF
3579	ENDIF
3580	#endif
3581
3582	ENDIF rk_step_1_check
3583
3584
3585	!**********************************************************
3586	!
3587	! end of RK predictor-corrector loop
3588	!
3589	!**********************************************************
3590
3591	END DO Runge_Kutta_loop
3592	#ifdef LMDZ1
3593	WRITE(message, *)' dyn_em: after runge_kutta_loop'
3594	CALL wrf_debug(200, message)
3595	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
3596	' u_tend: ', ru_tendf(im2,1,jm2)
3597	CALL wrf_debug(200, message)
3598	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
3599	'p sfc: ',p8w(im2,kms,jm2)
3600	CALL wrf_debug(200, message)
3601	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
3602	CALL wrf_debug(200, message)
3603	#endif
3604
3605	IF (config_flags%do_avgflx_em .EQ. 1) THEN
3606	! Reinitialize time-averaged fluxes if history output was written after the previous time step:
3607	CALL WRFU_ALARMGET(grid%alarms( HISTORY_ALARM ),prevringtime=temp_time)
3608	CALL domain_clock_get ( grid, current_time=CurrTime, &
3609	current_timestr=message2 )
3610	! use overloaded -, .LT. operator to check whether to initialize avgflx:
3611	! reinitialize after each history output (detect this here by comparing current time
3612	! against last history time and time step - this code follows what's done in adapt_timestep_em):
3613	WRITE ( message , FMT = '("solve_em: old_dt =",g15.6,", dt=",g15.6," on domain ",I3)' ) &
3614	& old_dt,grid%dt,grid%id
3615	CALL wrf_debug(200,message)
3616	old_dt=min(old_dt,grid%dt)
3617	num = INT(old_dt * precision)
3618	den = precision
3619	CALL WRFU_TimeIntervalSet(dtInterval, Sn=num, Sd=den)
3620	IF (CurrTime .lt. temp_time + dtInterval) THEN
3621	WRITE ( message , FMT = '("solve_em: initializing avgflx at time ",A," on domain ",I3)' ) &
3622	& TRIM(message2), grid%id
3623	CALL wrf_message(trim(message))
3624	grid%avgflx_count = 0
3625	!tile-loop for zero_avgflx
3626	!$OMP PARALLEL DO &
3627	!$OMP PRIVATE ( ij )
3628
3629	DO ij = 1 , grid%num_tiles
3630	CALL wrf_debug(200,'In solve_em, before zero_avgflx call')
3631	CALL zero_avgflx(grid%avgflx_rum,grid%avgflx_rvm,grid%avgflx_wwm, &
3632	& ids, ide, jds, jde, kds, kde, &
3633	& ims, ime, jms, jme, kms, kme, &
3634	& grid%i_start(ij), grid%i_end(ij), grid%j_start(ij), grid%j_end(ij), &
3635	& k_start , k_end, f_flux, &
3636	& grid%avgflx_cfu1,grid%avgflx_cfd1,grid%avgflx_dfu1, &
3637	& grid%avgflx_efu1,grid%avgflx_dfd1,grid%avgflx_efd1 )
3638	CALL wrf_debug(200,'In solve_em, after zero_avgflx call')
3639	ENDDO
3640	ENDIF
3641
3642	! Update avgflx quantities
3643	!tile-loop for upd_avgflx
3644	!$OMP PARALLEL DO &
3645	!$OMP PRIVATE ( ij )
3646
3647	DO ij = 1 , grid%num_tiles
3648	CALL wrf_debug(200,'In solve_em, before upd_avgflx call')
3649	CALL upd_avgflx(grid%avgflx_count,grid%avgflx_rum,grid%avgflx_rvm,grid%avgflx_wwm, &
3650	& grid%ru_m, grid%rv_m, grid%ww_m, &
3651	& ids, ide, jds, jde, kds, kde, &
3652	& ims, ime, jms, jme, kms, kme, &
3653	& grid%i_start(ij), grid%i_end(ij), grid%j_start(ij), grid%j_end(ij), &
3654	& k_start , k_end, f_flux, &
3655	& grid%cfu1,grid%cfd1,grid%dfu1,grid%efu1,grid%dfd1,grid%efd1, &
3656	& grid%avgflx_cfu1,grid%avgflx_cfd1,grid%avgflx_dfu1, &
3657	& grid%avgflx_efu1,grid%avgflx_dfd1,grid%avgflx_efd1 )
3658	CALL wrf_debug(200,'In solve_em, after upd_avgflx call')
3659
3660	ENDDO
3661	grid%avgflx_count = grid%avgflx_count + 1
3662	ENDIF
3663	!
3664	!$OMP PARALLEL DO &
3665	!$OMP PRIVATE ( ij )
3666	DO ij = 1 , grid%num_tiles
3667
3668	BENCH_START(advance_ppt_tim)
3669	CALL wrf_debug ( 200 , ' call advance_ppt' )
3670	CALL advance_ppt(grid%rthcuten,grid%rqvcuten,grid%rqccuten,grid%rqrcuten, &
3671	grid%rqicuten,grid%rqscuten, &
3672	grid%rainc,grid%raincv,grid%rainsh,grid%pratec,grid%pratesh, &
3673	grid%nca,grid%htop,grid%hbot,grid%cutop,grid%cubot, &
3674	grid%cuppt, grid%dt, config_flags, &
3675	ids,ide, jds,jde, kds,kde, &
3676	ims,ime, jms,jme, kms,kme, &
3677	grid%i_start(ij), grid%i_end(ij), &
3678	grid%j_start(ij), grid%j_end(ij), &
3679	k_start , k_end )
3680	BENCH_END(advance_ppt_tim)
3681
3682	ENDDO
3683	!$OMP END PARALLEL DO
3684
3685	!<DESCRIPTION>
3686	!<pre>
3687	! (5) time-split physics.
3688	!
3689	! Microphysics are the only time split physics in the WRF model
3690	! at this time. Split-physics begins with the calculation of
3691	! needed diagnostic quantities (pressure, temperature, etc.)
3692	! followed by a call to the microphysics driver,
3693	! and finishes with a clean-up, storing off of a diabatic tendency
3694	! from the moist physics, and a re-calulation of the diagnostic
3695	! quantities pressure and density.
3696	!</pre>
3697	!</DESCRIPTION>
3698
3699	IF( config_flags%specified .or. config_flags%nested ) THEN
3700	sz = grid%spec_zone
3701	ELSE
3702	sz = 0
3703	ENDIF
3704
3705	#ifdef LMDZ
3706	!!L. Fita, January 2105, LMD. We have mp_physics == 0, thus, WRF is not doing
3707	! this...
3708	PRINT *,' Lluis: num_3d_m: ',num_3d_m,' PARAM_FIRST_SCALAR: ', &
3709	PARAM_FIRST_SCALAR,' p_qv: ',p_qv,' p_qc: ',p_qc,' p_qr: ',p_qr, &
3710	' p_qi: ',p_qi,' p_qs: ',p_qs,' p_qg: ',p_qg, ' p_qh: ',p_qh
3711	IF (config_flags%mp_physics >= 0) then
3712	#else
3713	IF (config_flags%mp_physics /= 0) then
3714	#endif
3715	!$OMP PARALLEL DO &
3716	!$OMP PRIVATE ( ij, its, ite, jts, jte )
3717
3718	scalar_tile_loop_1a: DO ij = 1 , grid%num_tiles
3719
3720	IF ( config_flags%periodic_x ) THEN
3721	its = max(grid%i_start(ij),ids)
3722	ite = min(grid%i_end(ij),ide-1)
3723	ELSE
3724	its = max(grid%i_start(ij),ids+sz)
3725	ite = min(grid%i_end(ij),ide-1-sz)
3726	ENDIF
3727	jts = max(grid%j_start(ij),jds+sz)
3728	jte = min(grid%j_end(ij),jde-1-sz)
3729
3730	CALL wrf_debug ( 200 , ' call moist_physics_prep' )
3731	BENCH_START(moist_physics_prep_tim)
3732	CALL moist_physics_prep_em( grid%t_2, grid%t_1, t0, rho, &
3733	grid%al, grid%alb, grid%p, p8w, p0, grid%pb, &
3734	grid%ph_2, grid%phb, th_phy, pi_phy, p_phy, &
3735	grid%z, grid%z_at_w, dz8w, &
3736	dtm, grid%h_diabatic, &
3737	config_flags,grid%fnm, grid%fnp, &
3738	ids, ide, jds, jde, kds, kde, &
3739	ims, ime, jms, jme, kms, kme, &
3740	its, ite, jts, jte, &
3741	k_start , k_end )
3742	BENCH_END(moist_physics_prep_tim)
3743	END DO scalar_tile_loop_1a
3744	!$OMP END PARALLEL DO
3745
3746	CALL wrf_debug ( 200 , ' call microphysics_driver' )
3747
3748	grid%sr = 0.
3749	specified_bdy = config_flags%specified .OR. config_flags%nested
3750	channel_bdy = config_flags%specified .AND. config_flags%periodic_x
3751
3752	BENCH_START(micro_driver_tim)
3753
3754	#ifdef LMDZ
3755	! L. Fita, LMD February 2015. No moisture mircophysics!!
3756	! No microphyics_driver
3757	#else
3758	CALL microphysics_driver( &
3759	& DT=dtm ,DX=grid%dx ,DY=grid%dy &
3760	& ,DZ8W=dz8w ,F_ICE_PHY=grid%f_ice_phy &
3761	& ,ITIMESTEP=grid%itimestep ,LOWLYR=grid%lowlyr &
3762	& ,P8W=p8w ,P=p_phy ,PI_PHY=pi_phy &
3763	& ,RHO=rho ,SPEC_ZONE=grid%spec_zone &
3764	& ,SR=grid%sr ,TH=th_phy &
3765	& ,refl_10cm=grid%refl_10cm & ! hm, 9/22/09 for refl
3766	& ,WARM_RAIN=grid%warm_rain &
3767	& ,T8W=t8w &
3768	& ,CLDFRA=grid%cldfra, EXCH_H=grid%exch_h &
3769	& ,NSOURCE=grid%qndropsource &
3770	#ifdef WRF_CHEM
3771	& ,QLSINK=grid%qlsink,CLDFRA_OLD=grid%cldfra_old &
3772	& ,PRECR=grid%precr, PRECI=grid%preci, PRECS=grid%precs, PRECG=grid%precg &
3773	& ,CHEM_OPT=config_flags%chem_opt, PROGN=config_flags%progn &
3774	#endif
3775	& ,XLAND=grid%xland &
3776	& ,SPECIFIED=specified_bdy, CHANNEL_SWITCH=channel_bdy &
3777	& ,F_RAIN_PHY=grid%f_rain_phy &
3778	& ,F_RIMEF_PHY=grid%f_rimef_phy &
3779	& ,MP_PHYSICS=config_flags%mp_physics &
3780	& ,ID=grid%id &
3781	& ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
3782	& ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
3783	#ifdef RUN_ON_GPU
3784	& ,IPS=ips,IPE=ipe, JPS=jps,JPE=jpe, KPS=kps,KPE=kpe &
3785	#endif
3786	& ,I_START=grid%i_start,I_END=min(grid%i_end, ide-1) &
3787	& ,J_START=grid%j_start,J_END=min(grid%j_end, jde-1) &
3788	& ,KTS=k_start, KTE=min(k_end,kde-1) &
3789	& ,NUM_TILES=grid%num_tiles &
3790	& ,NAER=grid%naer &
3791	! Optional
3792	& , RAINNC=grid%rainnc, RAINNCV=grid%rainncv &
3793	& , SNOWNC=grid%snownc, SNOWNCV=grid%snowncv &
3794	& , GRAUPELNC=grid%graupelnc, GRAUPELNCV=grid%graupelncv & ! for milbrandt2mom
3795	& , HAILNC=grid%hailnc, HAILNCV=grid%hailncv &
3796	& , W=grid%w_2, Z=grid%z, HT=grid%ht &
3797	& , MP_RESTART_STATE=grid%mp_restart_state &
3798	& , TBPVS_STATE=grid%tbpvs_state & ! etampnew
3799	& , TBPVS0_STATE=grid%tbpvs0_state & ! etampnew
3800	& , QV_CURR=moist(ims,kms,jms,P_QV), F_QV=F_QV &
3801	& , QC_CURR=moist(ims,kms,jms,P_QC), F_QC=F_QC &
3802	& , QR_CURR=moist(ims,kms,jms,P_QR), F_QR=F_QR &
3803	& , QI_CURR=moist(ims,kms,jms,P_QI), F_QI=F_QI &
3804	& , QS_CURR=moist(ims,kms,jms,P_QS), F_QS=F_QS &
3805	& , QG_CURR=moist(ims,kms,jms,P_QG), F_QG=F_QG &
3806	& , QH_CURR=moist(ims,kms,jms,P_QH), F_QH=F_QH & ! for milbrandt2mom
3807	& , QNDROP_CURR=scalar(ims,kms,jms,P_QNDROP), F_QNDROP=F_QNDROP &
3808	& , QT_CURR=scalar(ims,kms,jms,P_QT), F_QT=F_QT &
3809	& , QNN_CURR=scalar(ims,kms,jms,P_QNN), F_QNN=F_QNN &
3810	& , QNI_CURR=scalar(ims,kms,jms,P_QNI), F_QNI=F_QNI &
3811	& , QNC_CURR=scalar(ims,kms,jms,P_QNC), F_QNC=F_QNC &
3812	& , QNR_CURR=scalar(ims,kms,jms,P_QNR), F_QNR=F_QNR &
3813	& , QNS_CURR=scalar(ims,kms,jms,P_QNS), F_QNS=F_QNS &
3814	& , QNG_CURR=scalar(ims,kms,jms,P_QNG), F_QNG=F_QNG &
3815	& , QNH_CURR=scalar(ims,kms,jms,P_QNH), F_QNH=F_QNH & ! for milbrandt2mom
3816	! & , QZR_CURR=scalar(ims,kms,jms,P_QZR), F_QZR=F_QZR & ! for milbrandt3mom
3817	! & , QZI_CURR=scalar(ims,kms,jms,P_QZI), F_QZI=F_QZI & ! "
3818	! & , QZS_CURR=scalar(ims,kms,jms,P_QZS), F_QZS=F_QZS & ! "
3819	! & , QZG_CURR=scalar(ims,kms,jms,P_QZG), F_QZG=F_QZG & ! "
3820	! & , QZH_CURR=scalar(ims,kms,jms,P_QZH), F_QZH=F_QZH & ! "
3821	& , qrcuten=grid%rqrcuten, qscuten=grid%rqscuten &
3822	& , qicuten=grid%rqicuten,mu=grid%mut &
3823	& , HAIL=config_flags%gsfcgce_hail & ! for gsfcgce
3824	& , ICE2=config_flags%gsfcgce_2ice & ! for gsfcgce
3825	! & , ccntype=config_flags%milbrandt_ccntype & ! for milbrandt (2mom)
3826	! YLIN
3827	! RI_CURR INPUT
3828	& , RI_CURR=grid%rimi &
3829	)
3830	#endif
3831	BENCH_END(micro_driver_tim)
3832
3833	#ifdef LMDZ
3834	! grid%h_diabatic = 0.
3835	#endif
3836
3837	#if 0
3838	BENCH_START(microswap_2)
3839	! for load balancing; communication to redistribute the points
3840	IF ( config_flags%mp_physics .EQ. ETAMPNEW ) THEN
3841	#include "SWAP_ETAMP_NEW.inc"
3842	ELSE IF ( config_flags%mp_physics .EQ. WSM3SCHEME ) THEN
3843	#include "SWAP_WSM3.inc"
3844	ENDIF
3845	BENCH_END(microswap_2)
3846	#endif
3847
3848	CALL wrf_debug ( 200 , ' call moist_physics_finish' )
3849	BENCH_START(moist_phys_end_tim)
3850
3851	!$OMP PARALLEL DO &
3852	!$OMP PRIVATE ( ij, its, ite, jts, jte, im, ii, jj, kk )
3853
3854	DO ij = 1 , grid%num_tiles
3855
3856	its = max(grid%i_start(ij),ids)
3857	ite = min(grid%i_end(ij),ide-1)
3858	jts = max(grid%j_start(ij),jds)
3859	jte = min(grid%j_end(ij),jde-1)
3860
3861	CALL microphysics_zero_outb ( &
3862	moist , num_moist , 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
3868	CALL microphysics_zero_outb ( &
3869	scalar , num_scalar , config_flags , &
3870	ids, ide, jds, jde, kds, kde, &
3871	ims, ime, jms, jme, kms, kme, &
3872	its, ite, jts, jte, &
3873	k_start , k_end )
3874
3875	CALL microphysics_zero_outb ( &
3876	chem , num_chem , config_flags , &
3877	ids, ide, jds, jde, kds, kde, &
3878	ims, ime, jms, jme, kms, kme, &
3879	its, ite, jts, jte, &
3880	k_start , k_end )
3881	CALL microphysics_zero_outb ( &
3882	tracer , num_tracer , config_flags , &
3883	ids, ide, jds, jde, kds, kde, &
3884	ims, ime, jms, jme, kms, kme, &
3885	its, ite, jts, jte, &
3886	k_start , k_end )
3887
3888	IF ( config_flags%periodic_x ) THEN
3889	its = max(grid%i_start(ij),ids)
3890	ite = min(grid%i_end(ij),ide-1)
3891	ELSE
3892	its = max(grid%i_start(ij),ids+sz)
3893	ite = min(grid%i_end(ij),ide-1-sz)
3894	ENDIF
3895	jts = max(grid%j_start(ij),jds+sz)
3896	jte = min(grid%j_end(ij),jde-1-sz)
3897
3898	CALL microphysics_zero_outa ( &
3899	moist , num_moist , 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	scalar , num_scalar , 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 microphysics_zero_outa ( &
3913	chem , num_chem , config_flags , &
3914	ids, ide, jds, jde, kds, kde, &
3915	ims, ime, jms, jme, kms, kme, &
3916	its, ite, jts, jte, &
3917	k_start , k_end )
3918
3919	CALL microphysics_zero_outa ( &
3920	tracer , num_tracer , config_flags , &
3921	ids, ide, jds, jde, kds, kde, &
3922	ims, ime, jms, jme, kms, kme, &
3923	its, ite, jts, jte, &
3924	k_start , k_end )
3925
3926	CALL moist_physics_finish_em( grid%t_2, grid%t_1, t0, grid%muts, th_phy, &
3927	grid%h_diabatic, dtm, config_flags, &
3928	#if ( WRF_DFI_RADAR == 1 )
3929	grid%dfi_tten_rad,grid%dfi_stage, &
3930	#endif
3931	ids, ide, jds, jde, kds, kde, &
3932	ims, ime, jms, jme, kms, kme, &
3933	its, ite, jts, jte, &
3934	k_start , k_end )
3935
3936	END DO
3937	!$OMP END PARALLEL DO
3938
3939	ENDIF ! microphysics test
3940
3941	!-----------------------------------------------------------
3942	! filter for moist variables post-microphysics and end of timestep
3943	!-----------------------------------------------------------
3944
3945	IF (config_flags%polar) THEN
3946	IF ( num_3d_m >= PARAM_FIRST_SCALAR ) THEN
3947	CALL wrf_debug ( 200 , ' call filter moist' )
3948	DO im = PARAM_FIRST_SCALAR, num_3d_m
3949	DO jj = jps, MIN(jpe,jde-1)
3950	DO kk = kps, MIN(kpe,kde-1)
3951	DO ii = ips, MIN(ipe,ide-1)
3952	moist(ii,kk,jj,im)=moist(ii,kk,jj,im)*(grid%mu_2(ii,jj)+grid%mub(ii,jj))
3953	ENDDO
3954	ENDDO
3955	ENDDO
3956
3957	CALL pxft ( grid=grid &
3958	,lineno=__LINE__ &
3959	,flag_uv = 0 &
3960	,flag_rurv = 0 &
3961	,flag_wph = 0 &
3962	,flag_ww = 0 &
3963	,flag_t = 0 &
3964	,flag_mu = 0 &
3965	,flag_mut = 0 &
3966	,flag_moist = im &
3967	,flag_chem = 0 &
3968	,flag_tracer = 0 &
3969	,flag_scalar = 0 &
3970	,positive_definite=.FALSE. &
3971	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
3972	,fft_filter_lat = config_flags%fft_filter_lat &
3973	,dclat = dclat &
3974	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
3975	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
3976	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
3977	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
3978	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
3979
3980	DO jj = jps, MIN(jpe,jde-1)
3981	DO kk = kps, MIN(kpe,kde-1)
3982	DO ii = ips, MIN(ipe,ide-1)
3983	moist(ii,kk,jj,im)=moist(ii,kk,jj,im)/(grid%mu_2(ii,jj)+grid%mub(ii,jj))
3984	ENDDO
3985	ENDDO
3986	ENDDO
3987	ENDDO
3988	ENDIF
3989	ENDIF
3990
3991	!-----------------------------------------------------------
3992	! end filter for moist variables post-microphysics and end of timestep
3993	!-----------------------------------------------------------
3994
3995	!$OMP PARALLEL DO &
3996	!$OMP PRIVATE ( ij, its, ite, jts, jte, im, ii, jj, kk )
3997	scalar_tile_loop_1ba: DO ij = 1 , grid%num_tiles
3998
3999	IF ( config_flags%periodic_x ) THEN
4000	its = max(grid%i_start(ij),ids)
4001	ite = min(grid%i_end(ij),ide-1)
4002	ELSE
4003	its = max(grid%i_start(ij),ids+sz)
4004	ite = min(grid%i_end(ij),ide-1-sz)
4005	ENDIF
4006	jts = max(grid%j_start(ij),jds+sz)
4007	jte = min(grid%j_end(ij),jde-1-sz)
4008
4009	CALL calc_p_rho_phi( moist, num_3d_m, &
4010	grid%al, grid%alb, grid%mu_2, grid%muts, &
4011	grid%ph_2, grid%p, grid%pb, grid%t_2, &
4012	p0, t0, grid%znu, grid%dnw, grid%rdnw, &
4013	grid%rdn, config_flags%non_hydrostatic, &
4014	ids, ide, jds, jde, kds, kde, &
4015	ims, ime, jms, jme, kms, kme, &
4016	its, ite, jts, jte, &
4017	k_start , k_end )
4018
4019	END DO scalar_tile_loop_1ba
4020	!$OMP END PARALLEL DO
4021	BENCH_END(moist_phys_end_tim)
4022
4023	IF (.not. config_flags%non_hydrostatic) THEN
4024	#ifdef DM_PARALLEL
4025	# include "HALO_EM_HYDRO_UV.inc"
4026	# include "PERIOD_EM_HYDRO_UV.inc"
4027	#endif
4028	!$OMP PARALLEL DO &
4029	!$OMP PRIVATE ( ij )
4030	DO ij = 1 , grid%num_tiles
4031	CALL diagnose_w( ph_tend, grid%ph_2, grid%ph_1, grid%w_2, grid%muts, dt_rk, &
4032	grid%u_2, grid%v_2, grid%ht, &
4033	grid%cf1, grid%cf2, grid%cf3, grid%rdx, grid%rdy, grid%msftx, grid%msfty, &
4034	ids, ide, jds, jde, kds, kde, &
4035	ims, ime, jms, jme, kms, kme, &
4036	grid%i_start(ij), grid%i_end(ij), &
4037	grid%j_start(ij), grid%j_end(ij), &
4038	k_start , k_end )
4039
4040	END DO
4041	!$OMP END PARALLEL DO
4042
4043	END IF
4044
4045	CALL wrf_debug ( 200 , ' call chem polar filter ' )
4046
4047	!-----------------------------------------------------------
4048	! filter for chem and scalar variables at end of timestep
4049	!-----------------------------------------------------------
4050
4051	IF (config_flags%polar) THEN
4052
4053	IF ( num_3d_c >= PARAM_FIRST_SCALAR ) then
4054	chem_filter_loop: DO im = PARAM_FIRST_SCALAR, num_3d_c
4055	DO jj = jps, MIN(jpe,jde-1)
4056	DO kk = kps, MIN(kpe,kde-1)
4057	DO ii = ips, MIN(ipe,ide-1)
4058	chem(ii,kk,jj,im)=chem(ii,kk,jj,im)*(grid%mu_2(ii,jj)+grid%mub(ii,jj))
4059	ENDDO
4060	ENDDO
4061	ENDDO
4062
4063	CALL pxft ( grid=grid &
4064	,lineno=__LINE__ &
4065	,flag_uv = 0 &
4066	,flag_rurv = 0 &
4067	,flag_wph = 0 &
4068	,flag_ww = 0 &
4069	,flag_t = 0 &
4070	,flag_mu = 0 &
4071	,flag_mut = 0 &
4072	,flag_moist = 0 &
4073	,flag_chem = im &
4074	,flag_tracer = 0 &
4075	,flag_scalar = 0 &
4076	,positive_definite=.FALSE. &
4077	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
4078	,fft_filter_lat = config_flags%fft_filter_lat &
4079	,dclat = dclat &
4080	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
4081	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
4082	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
4083	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
4084	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
4085
4086	DO jj = jps, MIN(jpe,jde-1)
4087	DO kk = kps, MIN(kpe,kde-1)
4088	DO ii = ips, MIN(ipe,ide-1)
4089	chem(ii,kk,jj,im)=chem(ii,kk,jj,im)/(grid%mu_2(ii,jj)+grid%mub(ii,jj))
4090	ENDDO
4091	ENDDO
4092	ENDDO
4093	ENDDO chem_filter_loop
4094	ENDIF
4095	IF ( num_tracer >= PARAM_FIRST_SCALAR ) then
4096	tracer_filter_loop: DO im = PARAM_FIRST_SCALAR, num_tracer
4097	DO jj = jps, MIN(jpe,jde-1)
4098	DO kk = kps, MIN(kpe,kde-1)
4099	DO ii = ips, MIN(ipe,ide-1)
4100	tracer(ii,kk,jj,im)=tracer(ii,kk,jj,im)*(grid%mu_2(ii,jj)+grid%mub(ii,jj))
4101	ENDDO
4102	ENDDO
4103	ENDDO
4104
4105	CALL pxft ( grid=grid &
4106	,lineno=__LINE__ &
4107	,flag_uv = 0 &
4108	,flag_rurv = 0 &
4109	,flag_wph = 0 &
4110	,flag_ww = 0 &
4111	,flag_t = 0 &
4112	,flag_mu = 0 &
4113	,flag_mut = 0 &
4114	,flag_moist = 0 &
4115	,flag_chem = 0 &
4116	,flag_tracer = im &
4117	,flag_scalar = 0 &
4118	,positive_definite=.FALSE. &
4119	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
4120	,fft_filter_lat = config_flags%fft_filter_lat &
4121	,dclat = dclat &
4122	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
4123	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
4124	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
4125	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
4126	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
4127
4128	DO jj = jps, MIN(jpe,jde-1)
4129	DO kk = kps, MIN(kpe,kde-1)
4130	DO ii = ips, MIN(ipe,ide-1)
4131	tracer(ii,kk,jj,im)=tracer(ii,kk,jj,im)/(grid%mu_2(ii,jj)+grid%mub(ii,jj))
4132	ENDDO
4133	ENDDO
4134	ENDDO
4135	ENDDO tracer_filter_loop
4136	ENDIF
4137
4138	IF ( num_3d_s >= PARAM_FIRST_SCALAR ) then
4139	scalar_filter_loop: DO im = PARAM_FIRST_SCALAR, num_3d_s
4140	DO jj = jps, MIN(jpe,jde-1)
4141	DO kk = kps, MIN(kpe,kde-1)
4142	DO ii = ips, MIN(ipe,ide-1)
4143	scalar(ii,kk,jj,im)=scalar(ii,kk,jj,im)*(grid%mu_2(ii,jj)+grid%mub(ii,jj))
4144	ENDDO
4145	ENDDO
4146	ENDDO
4147
4148	CALL pxft ( grid=grid &
4149	,lineno=__LINE__ &
4150	,flag_uv = 0 &
4151	,flag_rurv = 0 &
4152	,flag_wph = 0 &
4153	,flag_ww = 0 &
4154	,flag_t = 0 &
4155	,flag_mu = 0 &
4156	,flag_mut = 0 &
4157	,flag_moist = 0 &
4158	,flag_chem = 0 &
4159	,flag_tracer = 0 &
4160	,flag_scalar = im &
4161	,positive_definite=.FALSE. &
4162	,moist=moist,chem=chem,tracer=tracer,scalar=scalar &
4163	,fft_filter_lat = config_flags%fft_filter_lat &
4164	,dclat = dclat &
4165	,ids=ids,ide=ide,jds=jds,jde=jde,kds=kds,kde=kde &
4166	,ims=ims,ime=ime,jms=jms,jme=jme,kms=kms,kme=kme &
4167	,ips=ips,ipe=ipe,jps=jps,jpe=jpe,kps=kps,kpe=kpe &
4168	,imsx=imsx,imex=imex,jmsx=jmsx,jmex=jmex,kmsx=kmsx,kmex=kmex &
4169	,ipsx=ipsx,ipex=ipex,jpsx=jmsx,jpex=jpex,kpsx=kpsx,kpex=kpex )
4170
4171	DO jj = jps, MIN(jpe,jde-1)
4172	DO kk = kps, MIN(kpe,kde-1)
4173	DO ii = ips, MIN(ipe,ide-1)
4174	scalar(ii,kk,jj,im)=scalar(ii,kk,jj,im)/(grid%mu_2(ii,jj)+grid%mub(ii,jj))
4175	ENDDO
4176	ENDDO
4177	ENDDO
4178	ENDDO scalar_filter_loop
4179	ENDIF
4180	ENDIF
4181
4182	!-----------------------------------------------------------
4183	! end filter for chem and scalar variables at end of timestep
4184	!-----------------------------------------------------------
4185
4186	! We're finished except for boundary condition (and patch) update
4187
4188	! Boundary condition time (or communication time). At this time, we have
4189	! implemented periodic and symmetric physical boundary conditions.
4190
4191	! b.c. routine for data within patch.
4192
4193	! we need to do both time levels of
4194	! data because the time filter only works in the physical solution space.
4195
4196	! First, do patch communications for boundary conditions (periodicity)
4197
4198	!-----------------------------------------------------------
4199	! Stencils for patch communications (WCS, 29 June 2001)
4200	!
4201	! here's where we need a wide comm stencil - these are the
4202	! uncoupled variables so are used for high order calc in
4203	! advection and mixong routines.
4204	!
4205	! * * * * *
4206	! * * * * * * * * *
4207	! * + * * + * * * + * *
4208	! * * * * * * * * *
4209	! * * * * *
4210	!
4211	! grid%u_1 x
4212	! grid%u_2 x
4213	! grid%v_1 x
4214	! grid%v_2 x
4215	! grid%w_1 x
4216	! grid%w_2 x
4217	! grid%t_1 x
4218	! grid%t_2 x
4219	! grid%ph_1 x
4220	! grid%ph_2 x
4221	! grid%tke_1 x
4222	! grid%tke_2 x
4223	!
4224	! 2D variables
4225	! grid%mu_1 x
4226	! grid%mu_2 x
4227	!
4228	! 4D variables
4229	! moist x
4230	! chem x
4231	! scalar x
4232	!----------------------------------------------------------
4233
4234
4235	#ifdef DM_PARALLEL
4236	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
4237	# include "HALO_EM_D3_3.inc"
4238	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
4239	# include "HALO_EM_D3_5.inc"
4240	ELSE
4241	WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order
4242	CALL wrf_error_fatal(TRIM(wrf_err_message))
4243	ENDIF
4244	# include "PERIOD_BDY_EM_D3.inc"
4245	# include "PERIOD_BDY_EM_MOIST.inc"
4246	# include "PERIOD_BDY_EM_CHEM.inc"
4247	# include "PERIOD_BDY_EM_TRACER.inc"
4248	# include "PERIOD_BDY_EM_SCALAR.inc"
4249	#endif
4250
4251	! now set physical b.c on a patch
4252
4253	BENCH_START(bc_2d_tim)
4254	!$OMP PARALLEL DO &
4255	!$OMP PRIVATE ( ij )
4256	tile_bc_loop_2: DO ij = 1 , grid%num_tiles
4257
4258	CALL wrf_debug ( 200 , ' call set_phys_bc_dry_2' )
4259
4260	CALL set_phys_bc_dry_2( config_flags, &
4261	grid%u_1, grid%u_2, grid%v_1, grid%v_2, grid%w_1, grid%w_2, &
4262	grid%t_1, grid%t_2, grid%ph_1, grid%ph_2, grid%mu_1, grid%mu_2, &
4263	ids, ide, jds, jde, kds, kde, &
4264	ims, ime, jms, jme, kms, kme, &
4265	ips, ipe, jps, jpe, kps, kpe, &
4266	grid%i_start(ij), grid%i_end(ij), &
4267	grid%j_start(ij), grid%j_end(ij), &
4268	k_start , k_end )
4269
4270	CALL set_physical_bc3d( grid%tke_1, 'p', config_flags, &
4271	ids, ide, jds, jde, kds, kde, &
4272	ims, ime, jms, jme, kms, kme, &
4273	ips, ipe, jps, jpe, kps, kpe, &
4274	grid%i_start(ij), grid%i_end(ij), &
4275	grid%j_start(ij), grid%j_end(ij), &
4276	k_start , k_end-1 )
4277
4278	CALL set_physical_bc3d( grid%tke_2 , 'p', config_flags, &
4279	ids, ide, jds, jde, kds, kde, &
4280	ims, ime, jms, jme, kms, kme, &
4281	ips, ipe, jps, jpe, kps, kpe, &
4282	grid%i_start(ij), grid%i_end(ij), &
4283	grid%j_start(ij), grid%j_end(ij), &
4284	k_start , k_end )
4285
4286	moisture_loop_bdy_2 : DO im = PARAM_FIRST_SCALAR , num_3d_m
4287
4288	CALL set_physical_bc3d( moist(ims,kms,jms,im), 'p', &
4289	config_flags, &
4290	ids, ide, jds, jde, kds, kde, &
4291	ims, ime, jms, jme, kms, kme, &
4292	ips, ipe, jps, jpe, kps, kpe, &
4293	grid%i_start(ij), grid%i_end(ij), &
4294	grid%j_start(ij), grid%j_end(ij), &
4295	k_start , k_end )
4296
4297	END DO moisture_loop_bdy_2
4298
4299	chem_species_bdy_loop_2 : DO ic = PARAM_FIRST_SCALAR , num_3d_c
4300
4301	CALL set_physical_bc3d( chem(ims,kms,jms,ic) , 'p', config_flags, &
4302	ids, ide, jds, jde, kds, kde, &
4303	ims, ime, jms, jme, kms, kme, &
4304	ips, ipe, jps, jpe, kps, kpe, &
4305	grid%i_start(ij), grid%i_end(ij), &
4306	grid%j_start(ij), grid%j_end(ij), &
4307	k_start , k_end )
4308
4309	END DO chem_species_bdy_loop_2
4310
4311	tracer_species_bdy_loop_2 : DO ic = PARAM_FIRST_SCALAR , num_tracer
4312
4313	CALL set_physical_bc3d( tracer(ims,kms,jms,ic) , 'p', config_flags, &
4314	ids, ide, jds, jde, kds, kde, &
4315	ims, ime, jms, jme, kms, kme, &
4316	ips, ipe, jps, jpe, kps, kpe, &
4317	grid%i_start(ij), grid%i_end(ij), &
4318	grid%j_start(ij), grid%j_end(ij), &
4319	k_start , k_end )
4320
4321	END DO tracer_species_bdy_loop_2
4322
4323	scalar_species_bdy_loop_2 : DO is = PARAM_FIRST_SCALAR , num_3d_s
4324
4325	CALL set_physical_bc3d( scalar(ims,kms,jms,is) , 'p', config_flags, &
4326	ids, ide, jds, jde, kds, kde, &
4327	ims, ime, jms, jme, kms, kme, &
4328	ips, ipe, jps, jpe, kps, kpe, &
4329	grid%i_start(ij), grid%i_end(ij), &
4330	grid%j_start(ij), grid%j_end(ij), &
4331	k_start , k_end )
4332
4333	END DO scalar_species_bdy_loop_2
4334
4335	END DO tile_bc_loop_2
4336	!$OMP END PARALLEL DO
4337	BENCH_END(bc_2d_tim)
4338
4339	IF( config_flags%specified .or. config_flags%nested ) THEN
4340	grid%dtbc = grid%dtbc + grid%dt
4341	ENDIF
4342
4343	! reset surface w for consistency
4344
4345	#ifdef DM_PARALLEL
4346	# include "HALO_EM_C.inc"
4347	# include "PERIOD_BDY_EM_E.inc"
4348	#endif
4349
4350	CALL wrf_debug ( 10 , ' call set_w_surface' )
4351	fill_w_flag = .false.
4352
4353	!$OMP PARALLEL DO &
4354	!$OMP PRIVATE ( ij )
4355	DO ij = 1 , grid%num_tiles
4356	CALL set_w_surface( config_flags, grid%znw, fill_w_flag, &
4357	grid%w_2, grid%ht, grid%u_2, grid%v_2, &
4358	grid%cf1, grid%cf2, grid%cf3, grid%rdx, grid%rdy,&
4359	grid%msftx, grid%msfty, &
4360	ids, ide, jds, jde, kds, kde, &
4361	ims, ime, jms, jme, kms, kme, &
4362	grid%i_start(ij), grid%i_end(ij), &
4363	grid%j_start(ij), grid%j_end(ij), &
4364	k_start, k_end )
4365	! its, ite, jts, jte, k_start, min(k_end,kde-1), &
4366
4367	END DO
4368	!$OMP END PARALLEL DO
4369
4370	! calculate some model diagnostics.
4371
4372	#ifdef LMDZ1
4373	WRITE(message, *)' dyn_em: before diagnostics'
4374	CALL wrf_debug(200, message)
4375	WRITE(message, *)' t_tend: ',t_tendf(im2,km2,jm2), &
4376	' u_tend: ', ru_tendf(im2,1,jm2)
4377	CALL wrf_debug(200, message)
4378	WRITE(message,*)' psfc_tend: ',grid%dpsdt(im2,jm2), &
4379	'p sfc: ',p8w(im2,kms,jm2)
4380	CALL wrf_debug(200, message)
4381	WRITE(message,*)' p 1: ',grid%p(im2,kms,jm2), ' ph 1: ',grid%ph_2(im2,kms,jm2)
4382	CALL wrf_debug(200, message)
4383	#endif
4384	CALL wrf_debug ( 200 , ' call diagnostic_driver' )
4385
4386	CALL diagnostic_output_calc( &
4387	& DPSDT=grid%dpsdt ,DMUDT=grid%dmudt &
4388	& ,P8W=p8w ,PK1M=grid%pk1m &
4389	& ,MU_2=grid%mu_2 ,MU_2M=grid%mu_2m &
4390	& ,U=grid%u_2 ,V=grid%v_2 &
4391	& ,RAINCV=grid%raincv ,RAINNCV=grid%rainncv &
4392	& ,RAINC=grid%rainc ,RAINNC=grid%rainnc &
4393	& ,I_RAINC=grid%i_rainc ,I_RAINNC=grid%i_rainnc &
4394	& ,HFX=grid%hfx ,SFCEVP=grid%sfcevp ,LH=grid%lh &
4395	& ,DT=grid%dt ,SBW=config_flags%spec_bdy_width &
4396	& ,XTIME=grid%xtime ,T2=grid%t2 &
4397	& ,ACSWUPT=grid%acswupt ,ACSWUPTC=grid%acswuptc &
4398	& ,ACSWDNT=grid%acswdnt ,ACSWDNTC=grid%acswdntc &
4399	& ,ACSWUPB=grid%acswupb ,ACSWUPBC=grid%acswupbc &
4400	& ,ACSWDNB=grid%acswdnb ,ACSWDNBC=grid%acswdnbc &
4401	& ,ACLWUPT=grid%aclwupt ,ACLWUPTC=grid%aclwuptc &
4402	& ,ACLWDNT=grid%aclwdnt ,ACLWDNTC=grid%aclwdntc &
4403	& ,ACLWUPB=grid%aclwupb ,ACLWUPBC=grid%aclwupbc &
4404	& ,ACLWDNB=grid%aclwdnb ,ACLWDNBC=grid%aclwdnbc &
4405	& ,I_ACSWUPT=grid%i_acswupt ,I_ACSWUPTC=grid%i_acswuptc &
4406	& ,I_ACSWDNT=grid%i_acswdnt ,I_ACSWDNTC=grid%i_acswdntc &
4407	& ,I_ACSWUPB=grid%i_acswupb ,I_ACSWUPBC=grid%i_acswupbc &
4408	& ,I_ACSWDNB=grid%i_acswdnb ,I_ACSWDNBC=grid%i_acswdnbc &
4409	& ,I_ACLWUPT=grid%i_aclwupt ,I_ACLWUPTC=grid%i_aclwuptc &
4410	& ,I_ACLWDNT=grid%i_aclwdnt ,I_ACLWDNTC=grid%i_aclwdntc &
4411	& ,I_ACLWUPB=grid%i_aclwupb ,I_ACLWUPBC=grid%i_aclwupbc &
4412	& ,I_ACLWDNB=grid%i_aclwdnb ,I_ACLWDNBC=grid%i_aclwdnbc &
4413	! Selection flag
4414	& ,DIAG_PRINT=config_flags%diag_print &
4415	& ,BUCKET_MM=config_flags%bucket_mm &
4416	& ,BUCKET_J =config_flags%bucket_J &
4417	& ,SNOWNCV=grid%snowncv, SNOW_ACC_NC=grid%snow_acc_nc &
4418	& ,PREC_ACC_C=grid%prec_acc_c &
4419	& ,PREC_ACC_NC=grid%prec_acc_nc &
4420	& ,PREC_ACC_DT=config_flags%prec_acc_dt &
4421	& ,CURR_SECS=curr_secs &
4422	! Dimension arguments
4423	& ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
4424	& ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
4425	& ,IPS=ips,IPE=ipe, JPS=jps,JPE=jpe, KPS=kps,KPE=kpe &
4426	& ,I_START=grid%i_start,I_END=min(grid%i_end, ide-1) &
4427	& ,J_START=grid%j_start,J_END=min(grid%j_end, jde-1) &
4428	& ,KTS=k_start, KTE=min(k_end,kde-1) &
4429	& ,NUM_TILES=grid%num_tiles &
4430	& )
4431
4432	#ifdef DM_PARALLEL
4433	!-----------------------------------------------------------------------
4434	! see above
4435	!--------------------------------------------------------------
4436	CALL wrf_debug ( 200 , ' call HALO_RK_E' )
4437	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
4438	# include "HALO_EM_E_3.inc"
4439	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
4440	# include "HALO_EM_E_5.inc"
4441	ELSE
4442	WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order
4443	CALL wrf_error_fatal(TRIM(wrf_err_message))
4444	ENDIF
4445	#endif
4446
4447	#ifdef DM_PARALLEL
4448	IF ( num_moist >= PARAM_FIRST_SCALAR ) THEN
4449	!-----------------------------------------------------------------------
4450	! see above
4451	!--------------------------------------------------------------
4452	CALL wrf_debug ( 200 , ' call HALO_RK_MOIST' )
4453	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
4454	# include "HALO_EM_MOIST_E_3.inc"
4455	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
4456	# include "HALO_EM_MOIST_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_chem >= PARAM_FIRST_SCALAR ) THEN
4463	!-----------------------------------------------------------------------
4464	! see above
4465	!--------------------------------------------------------------
4466	CALL wrf_debug ( 200 , ' call HALO_RK_CHEM' )
4467	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
4468	# include "HALO_EM_CHEM_E_3.inc"
4469	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
4470	# include "HALO_EM_CHEM_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_tracer >= PARAM_FIRST_SCALAR ) THEN
4477	!-----------------------------------------------------------------------
4478	! see above
4479	!--------------------------------------------------------------
4480	CALL wrf_debug ( 200 , ' call HALO_RK_TRACER' )
4481	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
4482	# include "HALO_EM_TRACER_E_3.inc"
4483	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
4484	# include "HALO_EM_TRACER_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	IF ( num_scalar >= PARAM_FIRST_SCALAR ) THEN
4491	!-----------------------------------------------------------------------
4492	! see above
4493	!--------------------------------------------------------------
4494	CALL wrf_debug ( 200 , ' call HALO_RK_SCALAR' )
4495	IF ( config_flags%h_mom_adv_order <= 4 ) THEN
4496	# include "HALO_EM_SCALAR_E_3.inc"
4497	ELSE IF ( config_flags%h_mom_adv_order <= 6 ) THEN
4498	# include "HALO_EM_SCALAR_E_5.inc"
4499	ELSE
4500	WRITE(wrf_err_message,*)'solve_em: invalid h_mom_adv_order = ',config_flags%h_mom_adv_order
4501	CALL wrf_error_fatal(TRIM(wrf_err_message))
4502	ENDIF
4503	ENDIF
4504	#endif
4505
4506	! Max values of CFL for adaptive time step scheme
4507
4508	DEALLOCATE(max_vert_cfl_tmp)
4509	DEALLOCATE(max_horiz_cfl_tmp)
4510
4511	CALL wrf_debug ( 200 , ' call end of solve_em' )
4512
4513	! Finish timers if compiled with -DBENCH.
4514	#include <bench_solve_em_end.h>
4515
4516	RETURN
4517
4518	END SUBROUTINE solve_em

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

Download in other formats:

Original Format