Home

Context Navigation

source: palm/trunk/SOURCE/time_integration.f90 @ 2968

Last change on this file since 2968 was 2941, checked in by kanani, 6 years ago
Fixes for spinup and sky-view-factor (SVF) treatment
Property svn:keywords set to `Id`
File size: 51.4 KB

Line
1	!> @file time_integration.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2018 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: time_integration.f90 2941 2018-04-03 11:54:58Z suehring $
27	! Deduct spinup_time from RUN_CONTROL output of main 3d run
28	! (use time_since_reference_point instead of simulated_time)
29	!
30	! 2938 2018-03-27 15:52:42Z suehring
31	! Nesting of dissipation rate in case of RANS mode and TKE-e closure is applied
32	!
33	! 2936 2018-03-27 14:49:27Z suehring
34	! Little formatting adjustment.
35	!
36	! 2817 2018-02-19 16:32:21Z knoop
37	! Preliminary gust module interface implemented
38	!
39	! 2801 2018-02-14 16:01:55Z thiele
40	! Changed lpm from subroutine to module.
41	! Introduce particle transfer in nested models.
42	!
43	! 2776 2018-01-31 10:44:42Z Giersch
44	! Variable use_synthetic_turbulence_generator has been abbreviated
45	!
46	! 2773 2018-01-30 14:12:54Z suehring
47	! - Nesting for chemical species
48	!
49	! 2766 2018-01-22 17:17:47Z kanani
50	! Removed preprocessor directive __chem
51	!
52	! 2718 2018-01-02 08:49:38Z maronga
53	! Corrected "Former revisions" section
54	!
55	! 2696 2017-12-14 17:12:51Z kanani
56	! - Change in file header (GPL part)
57	! - Implementation of uv exposure model (FK)
58	! - Moved vnest_boundary_conds_khkm from tcm_diffusivities to here (TG)
59	! - renamed diffusivities to tcm_diffusivities (TG)
60	! - implement prognostic equation for diss (TG)
61	! - Moved/commented CALL to chem_emissions (FK)
62	! - Added CALL to chem_emissions (FK)
63	! - Implementation of chemistry module (FK)
64	! - Calls for setting boundary conditions in USM and LSM (MS)
65	! - Large-scale forcing with larger-scale models implemented (MS)
66	! - Rename usm_radiation into radiation_interactions; merge with branch
67	! radiation (MS)
68	! - added call for usm_green_heat_model for green building surfaces (RvT)
69	! - added call for usm_temperature_near_surface for use in indoor model (RvT)
70	!
71	! 2617 2017-11-16 12:47:24Z suehring
72	! Bugfix, assure that the reference state does not become zero.
73	!
74	! 2563 2017-10-19 15:36:10Z Giersch
75	! Variable wind_turbine moved to module control_parameters
76	!
77	! 2365 2017-08-21 14:59:59Z kanani
78	! Vertical grid nesting implemented (SadiqHuq)
79	!
80	! 2320 2017-07-21 12:47:43Z suehring
81	! Set bottom boundary conditions after nesting interpolation and anterpolation
82	!
83	! 2299 2017-06-29 10:14:38Z maronga
84	! Call of soil model adjusted
85	!
86	! 2292 2017-06-20 09:51:42Z schwenkel
87	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
88	! includes two more prognostic equations for cloud drop concentration (nc)
89	! and cloud water content (qc).
90	!
91	! 2271 2017-06-09 12:34:55Z sward
92	! Start timestep message changed
93	!
94	! 2259 2017-06-08 09:09:11Z gronemeier
95	! Implemented synthetic turbulence generator
96	!
97	! 2233 2017-05-30 18:08:54Z suehring
98	!
99	! 2232 2017-05-30 17:47:52Z suehring
100	! Adjustments to new topography and surface concept
101	! Modify passed parameters for disturb_field
102	!
103	! 2178 2017-03-17 11:07:39Z hellstea
104	! Setting perturbations at all times near inflow boundary is removed
105	! in case of nested boundaries
106	!
107	! 2174 2017-03-13 08:18:57Z maronga
108	! Added support for nesting with cloud microphysics
109	!
110	! 2118 2017-01-17 16:38:49Z raasch
111	! OpenACC directives and related code removed
112	!
113	! 2050 2016-11-08 15:00:55Z gronemeier
114	! Implement turbulent outflow condition
115	!
116	! 2031 2016-10-21 15:11:58Z knoop
117	! renamed variable rho to rho_ocean
118	!
119	! 2011 2016-09-19 17:29:57Z kanani
120	! Flag urban_surface is now defined in module control_parameters,
121	! removed commented CALLs of global_min_max.
122	!
123	! 2007 2016-08-24 15:47:17Z kanani
124	! Added CALLs for new urban surface model
125	!
126	! 2000 2016-08-20 18:09:15Z knoop
127	! Forced header and separation lines into 80 columns
128	!
129	! 1976 2016-07-27 13:28:04Z maronga
130	! Simplified calls to radiation model
131	!
132	! 1960 2016-07-12 16:34:24Z suehring
133	! Separate humidity and passive scalar
134	!
135	! 1957 2016-07-07 10:43:48Z suehring
136	! flight module added
137	!
138	! 1919 2016-05-27 14:51:23Z raasch
139	! Initial version of purely vertical nesting introduced.
140	!
141	! 1918 2016-05-27 14:35:57Z raasch
142	! determination of time step moved to the end of the time step loop,
143	! the first time step is now always calculated before the time step loop (i.e.
144	! also in case of restart runs)
145	!
146	! 1914 2016-05-26 14:44:07Z witha
147	! Added call for wind turbine model
148	!
149	! 1878 2016-04-19 12:30:36Z hellstea
150	! Synchronization for nested runs rewritten
151	!
152	! 1853 2016-04-11 09:00:35Z maronga
153	! Adjusted for use with radiation_scheme = constant
154	!
155	! 1849 2016-04-08 11:33:18Z hoffmann
156	! Adapted for modularization of microphysics
157	!
158	! 1833 2016-04-07 14:23:03Z raasch
159	! spectrum renamed spectra_mod, spectra related variables moved to spectra_mod
160	!
161	! 1831 2016-04-07 13:15:51Z hoffmann
162	! turbulence renamed collision_turbulence
163	!
164	! 1822 2016-04-07 07:49:42Z hoffmann
165	! icloud_scheme replaced by microphysics_*
166	!
167	! 1808 2016-04-05 19:44:00Z raasch
168	! output message in case unscheduled radiation calls removed
169	!
170	! 1797 2016-03-21 16:50:28Z raasch
171	! introduction of different datatransfer modes
172	!
173	! 1791 2016-03-11 10:41:25Z raasch
174	! call of pmci_update_new removed
175	!
176	! 1786 2016-03-08 05:49:27Z raasch
177	! +module spectrum
178	!
179	! 1783 2016-03-06 18:36:17Z raasch
180	! switch back of netcdf data format for mask output moved to the mask output
181	! routine
182	!
183	! 1781 2016-03-03 15:12:23Z raasch
184	! some pmc calls removed at the beginning (before timeloop),
185	! pmc initialization moved to the main program
186	!
187	! 1764 2016-02-28 12:45:19Z raasch
188	! PMC_ACTIVE flags removed,
189	! bugfix: nest synchronization after first call of timestep
190	!
191	! 1762 2016-02-25 12:31:13Z hellstea
192	! Introduction of nested domain feature
193	!
194	! 1736 2015-12-04 08:56:33Z raasch
195	! no perturbations added to total domain if energy limit has been set zero
196	!
197	! 1691 2015-10-26 16:17:44Z maronga
198	! Added option for spin-ups without land surface and radiation models. Moved calls
199	! for radiation and lan surface schemes.
200	!
201	! 1682 2015-10-07 23:56:08Z knoop
202	! Code annotations made doxygen readable
203	!
204	! 1671 2015-09-25 03:29:37Z raasch
205	! bugfix: ghostpoint exchange for array diss in case that sgs velocities are used
206	! for particles
207	!
208	! 1585 2015-04-30 07:05:52Z maronga
209	! Moved call of radiation scheme. Added support for RRTM
210	!
211	! 1551 2015-03-03 14:18:16Z maronga
212	! Added interface for different radiation schemes.
213	!
214	! 1496 2014-12-02 17:25:50Z maronga
215	! Added calls for the land surface model and radiation scheme
216	!
217	! 1402 2014-05-09 14:25:13Z raasch
218	! location messages modified
219	!
220	! 1384 2014-05-02 14:31:06Z raasch
221	! location messages added
222	!
223	! 1380 2014-04-28 12:40:45Z heinze
224	! CALL of nudge_ref added
225	! bc_pt_t_val and bc_q_t_val are updated in case nudging is used
226	!
227	! 1365 2014-04-22 15:03:56Z boeske
228	! Reset sums_ls_l to zero at each timestep
229	! +sums_ls_l
230	! Calculation of reference state (previously in subroutine calc_mean_profile)
231
232	! 1342 2014-03-26 17:04:47Z kanani
233	! REAL constants defined as wp-kind
234	!
235	! 1320 2014-03-20 08:40:49Z raasch
236	! ONLY-attribute added to USE-statements,
237	! kind-parameters added to all INTEGER and REAL declaration statements,
238	! kinds are defined in new module kinds,
239	! old module precision_kind is removed,
240	! revision history before 2012 removed,
241	! comment fields (!:) to be used for variable explanations added to
242	! all variable declaration statements
243	! 1318 2014-03-17 13:35:16Z raasch
244	! module interfaces removed
245	!
246	! 1308 2014-03-13 14:58:42Z fricke
247	! +netcdf_data_format_save
248	! For masked data, parallel netcdf output is not tested so far, hence
249	! netcdf_data_format is switched back to non-paralell output.
250	!
251	! 1276 2014-01-15 13:40:41Z heinze
252	! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
253	!
254	! 1257 2013-11-08 15:18:40Z raasch
255	! acc-update-host directive for timestep removed
256	!
257	! 1241 2013-10-30 11:36:58Z heinze
258	! Generalize calc_mean_profile for wider use
259	! Determine shf and qsws in dependence on data from LSF_DATA
260	! Determine ug and vg in dependence on data from LSF_DATA
261	! 1221 2013-09-10 08:59:13Z raasch
262	! host update of arrays before timestep is called
263	!
264	! 1179 2013-06-14 05:57:58Z raasch
265	! mean profiles for reference state are only calculated if required,
266	! small bugfix for background communication
267	!
268	! 1171 2013-05-30 11:27:45Z raasch
269	! split of prognostic_equations deactivated (comment lines), for the time being
270	!
271	! 1128 2013-04-12 06:19:32Z raasch
272	! asynchronous transfer of ghost point data realized for acc-optimized version:
273	! prognostic_equations are first called two times for those points required for
274	! the left-right and north-south exchange, respectively, and then for the
275	! remaining points,
276	! those parts requiring global communication moved from prognostic_equations to
277	! here
278	!
279	! 1115 2013-03-26 18:16:16Z hoffmann
280	! calculation of qr and nr is restricted to precipitation
281	!
282	! 1113 2013-03-10 02:48:14Z raasch
283	! GPU-porting of boundary conditions,
284	! openACC directives updated
285	! formal parameter removed from routine boundary_conds
286	!
287	! 1111 2013-03-08 23:54:10Z raasch
288	! +internal timestep counter for cpu statistics added,
289	! openACC directives updated
290	!
291	! 1092 2013-02-02 11:24:22Z raasch
292	! unused variables removed
293	!
294	! 1065 2012-11-22 17:42:36Z hoffmann
295	! exchange of diss (dissipation rate) in case of turbulence = .TRUE. added
296	!
297	! 1053 2012-11-13 17:11:03Z hoffmann
298	! exchange of ghost points for nr, qr added
299	!
300	! 1036 2012-10-22 13:43:42Z raasch
301	! code put under GPL (PALM 3.9)
302	!
303	! 1019 2012-09-28 06:46:45Z raasch
304	! non-optimized version of prognostic_equations removed
305	!
306	! 1015 2012-09-27 09:23:24Z raasch
307	! +call of prognostic_equations_acc
308	!
309	! 1001 2012-09-13 14:08:46Z raasch
310	! all actions concerning leapfrog- and upstream-spline-scheme removed
311	!
312	! 849 2012-03-15 10:35:09Z raasch
313	! advec_particles renamed lpm, first_call_advec_particles renamed first_call_lpm
314	!
315	! 825 2012-02-19 03:03:44Z raasch
316	! wang_collision_kernel renamed wang_kernel
317	!
318	! Revision 1.1 1997/08/11 06:19:04 raasch
319	! Initial revision
320	!
321	!
322	! Description:
323	! ------------
324	!> Integration in time of the model equations, statistical analysis and graphic
325	!> output
326	!------------------------------------------------------------------------------!
327	SUBROUTINE time_integration
328
329
330	USE advec_ws, &
331	ONLY: ws_statistics
332
333	USE arrays_3d, &
334	ONLY: diss, diss_p, dzu, e, e_p, nc, nc_p, nr, nr_p, prho, pt, pt_p, pt_init, &
335	q_init, q, qc, qc_p, ql, ql_c, ql_v, ql_vp, qr, qr_p, q_p, &
336	ref_state, rho_ocean, s, s_p, sa_p, tend, u, u_p, v, vpt, &
337	v_p, w, w_p
338
339	USE calc_mean_profile_mod, &
340	ONLY: calc_mean_profile
341
342	USE chemistry_model_mod, &
343	ONLY: chem_emissions, chem_species
344
345	USE chem_modules, &
346	ONLY: nspec
347
348	USE control_parameters, &
349	ONLY: advected_distance_x, advected_distance_y, air_chemistry, &
350	average_count_3d, averaging_interval, averaging_interval_pr, &
351	bc_lr_cyc, bc_ns_cyc, bc_pt_t_val, &
352	bc_q_t_val, call_psolver_at_all_substeps, cloud_droplets, &
353	cloud_physics, constant_flux_layer, constant_heatflux, &
354	create_disturbances, dopr_n, constant_diffusion, coupling_mode, &
355	coupling_start_time, current_timestep_number, &
356	disturbance_created, disturbance_energy_limit, dist_range, &
357	do_sum, dt_3d, dt_averaging_input, dt_averaging_input_pr, &
358	dt_coupling, dt_data_output_av, dt_disturb, dt_do2d_xy, &
359	dt_do2d_xz, dt_do2d_yz, dt_do3d, dt_domask,dt_dopts, dt_dopr, &
360	dt_dopr_listing, dt_dots, dt_dvrp, dt_run_control, end_time, &
361	first_call_lpm, forcing, galilei_transformation, humidity, &
362	intermediate_timestep_count, intermediate_timestep_count_max, &
363	land_surface, large_scale_forcing, &
364	loop_optimization, lsf_surf, lsf_vert, masks, &
365	microphysics_morrison, microphysics_seifert, mid, nest_domain, &
366	neutral, nr_timesteps_this_run, nudging, &
367	ocean, passive_scalar, prho_reference, pt_reference, &
368	pt_slope_offset, random_heatflux, rans_mode, &
369	rans_tke_e, run_coupled, simulated_time, simulated_time_chr, &
370	skip_time_do2d_xy, skip_time_do2d_xz, skip_time_do2d_yz, &
371	skip_time_do3d, skip_time_domask, skip_time_dopr, &
372	skip_time_data_output_av, sloping_surface, &
373	stop_dt, terminate_coupled, terminate_run, timestep_scheme, &
374	time_coupling, time_do2d_xy, time_do2d_xz, time_do2d_yz, &
375	time_do3d, time_domask, time_dopr, time_dopr_av, &
376	time_dopr_listing, time_dopts, time_dosp, time_dosp_av, &
377	time_dots, time_do_av, time_do_sla, time_disturb, time_dvrp, &
378	time_run_control, time_since_reference_point, &
379	turbulent_inflow, turbulent_outflow, urban_surface, &
380	use_initial_profile_as_reference, &
381	use_single_reference_value, uv_exposure, u_gtrans, v_gtrans, &
382	virtual_flight, wind_turbine, ws_scheme_mom, ws_scheme_sca
383
384	USE cpulog, &
385	ONLY: cpu_log, log_point, log_point_s
386
387	USE flight_mod, &
388	ONLY: flight_measurement
389
390	USE gust_mod, &
391	ONLY: gust_actions, gust_module_enabled
392
393	USE indices, &
394	ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
395
396	USE interaction_droplets_ptq_mod, &
397	ONLY: interaction_droplets_ptq
398
399	USE interfaces
400
401	USE kinds
402
403	USE land_surface_model_mod, &
404	ONLY: lsm_boundary_condition, lsm_energy_balance, lsm_soil_model, &
405	skip_time_do_lsm
406
407	USE lsf_nudging_mod, &
408	ONLY: calc_tnudge, ls_forcing_surf, ls_forcing_vert, nudge_ref, &
409	forcing_bc, forcing_bc_mass_conservation
410
411	USE netcdf_data_input_mod, &
412	ONLY: force, netcdf_data_input_lsf
413
414	USE microphysics_mod, &
415	ONLY: collision_turbulence
416
417	USE particle_attributes, &
418	ONLY: particle_advection, particle_advection_start, &
419	use_sgs_for_particles, wang_kernel
420
421	USE pegrid
422
423	USE pmc_interface, &
424	ONLY: nested_run, nesting_mode, pmci_boundary_conds, pmci_datatrans, &
425	pmci_ensure_nest_mass_conservation, pmci_synchronize
426
427	USE progress_bar, &
428	ONLY: finish_progress_bar, output_progress_bar
429
430	USE prognostic_equations_mod, &
431	ONLY: prognostic_equations_cache, prognostic_equations_vector
432
433	USE radiation_model_mod, &
434	ONLY: dt_radiation, force_radiation_call, radiation, radiation_control,&
435	radiation_interaction, radiation_interactions, &
436	skip_time_do_radiation, time_radiation
437
438	USE spectra_mod, &
439	ONLY: average_count_sp, averaging_interval_sp, calc_spectra, dt_dosp, &
440	skip_time_dosp
441
442	USE statistics, &
443	ONLY: flow_statistics_called, hom, pr_palm, sums_ls_l, u_max, &
444	u_max_ijk, v_max, v_max_ijk, w_max, w_max_ijk
445
446	USE surface_layer_fluxes_mod, &
447	ONLY: surface_layer_fluxes
448
449	USE surface_mod, &
450	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
451
452	USE turbulence_closure_mod, &
453	ONLY: tcm_diffusivities, production_e_init
454
455	USE urban_surface_mod, &
456	ONLY: usm_boundary_condition, usm_material_heat_model, &
457	usm_material_model, &
458	usm_surface_energy_balance, usm_green_heat_model, &
459	usm_temperature_near_surface
460
461	USE synthetic_turbulence_generator_mod, &
462	ONLY: stg_main, use_syn_turb_gen
463
464	USE user_actions_mod, &
465	ONLY: user_actions
466
467	USE uv_exposure_model_mod, &
468	ONLY: uvem_calc_exposure
469
470	USE wind_turbine_model_mod, &
471	ONLY: wtm_forces
472
473	USE lpm_mod, &
474	ONLY: lpm
475
476	USE vertical_nesting_mod, &
477	ONLY: vnested, vnest_anterpolate, vnest_anterpolate_e, &
478	vnest_boundary_conds, vnest_boundary_conds_khkm, &
479	vnest_deallocate, vnest_init, vnest_init_fine, &
480	vnest_start_time
481
482	IMPLICIT NONE
483
484	CHARACTER (LEN=9) :: time_to_string !<
485	INTEGER :: n
486	INTEGER :: lsp
487
488	REAL(wp) :: dt_3d_old !< temporary storage of timestep to be used for
489	!< steering of run control output interval
490	!
491	!-- At beginning determine the first time step
492	CALL timestep
493	!
494	!-- Synchronize the timestep in case of nested run.
495	IF ( nested_run ) THEN
496	!
497	!-- Synchronization by unifying the time step.
498	!-- Global minimum of all time-steps is used for all.
499	CALL pmci_synchronize
500	ENDIF
501
502	!
503	!-- Determine and print out the run control quantities before the first time
504	!-- step of this run. For the initial run, some statistics (e.g. divergence)
505	!-- need to be determined first.
506	IF ( simulated_time == 0.0_wp ) CALL flow_statistics
507	CALL run_control
508	!
509	!-- Data exchange between coupled models in case that a call has been omitted
510	!-- at the end of the previous run of a job chain.
511	IF ( coupling_mode /= 'uncoupled' .AND. run_coupled .AND. .NOT. vnested) THEN
512	!
513	!-- In case of model termination initiated by the local model the coupler
514	!-- must not be called because this would again cause an MPI hang.
515	DO WHILE ( time_coupling >= dt_coupling .AND. terminate_coupled == 0 )
516	CALL surface_coupler
517	time_coupling = time_coupling - dt_coupling
518	ENDDO
519	IF (time_coupling == 0.0_wp .AND. &
520	time_since_reference_point < dt_coupling ) &
521	THEN
522	time_coupling = time_since_reference_point
523	ENDIF
524	ENDIF
525
526	#if defined( __dvrp_graphics )
527	!
528	!-- Time measurement with dvrp software
529	CALL DVRP_LOG_EVENT( 2, current_timestep_number )
530	#endif
531
532	CALL location_message( 'starting timestep-sequence', .TRUE. )
533	!
534	!-- Start of the time loop
535	DO WHILE ( simulated_time < end_time .AND. .NOT. stop_dt .AND. &
536	.NOT. terminate_run )
537
538	CALL cpu_log( log_point_s(10), 'timesteps', 'start' )
539	!
540	!-- Vertical nesting: initialize fine grid
541	IF ( vnested ) THEN
542	IF ( .NOT. vnest_init .AND. simulated_time >= vnest_start_time ) THEN
543	CALL cpu_log( log_point(80), 'vnest_init', 'start' )
544	CALL vnest_init_fine
545	vnest_init = .TRUE.
546	CALL cpu_log( log_point(80), 'vnest_init', 'stop' )
547	ENDIF
548	ENDIF
549	!
550	!-- Determine ug, vg and w_subs in dependence on data from external file
551	!-- LSF_DATA
552	IF ( large_scale_forcing .AND. lsf_vert ) THEN
553	CALL ls_forcing_vert ( simulated_time )
554	sums_ls_l = 0.0_wp
555	ENDIF
556
557	!
558	!-- Set pt_init and q_init to the current profiles taken from
559	!-- NUDGING_DATA
560	IF ( nudging ) THEN
561	CALL nudge_ref ( simulated_time )
562	!
563	!-- Store temperature gradient at the top boundary for possible Neumann
564	!-- boundary condition
565	bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
566	bc_q_t_val = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
567	ENDIF
568	!
569	!-- If forcing by larger-scale models is applied, check if new data
570	!-- at domain boundaries need to be read.
571	IF ( forcing ) THEN
572	IF ( force%time(force%tind_p) <= simulated_time ) &
573	CALL netcdf_data_input_lsf
574	ENDIF
575
576	!
577	!-- Execute the gust module actions
578	IF ( gust_module_enabled ) THEN
579	CALL gust_actions( 'before_timestep' )
580	ENDIF
581
582	!
583	!-- Execute the user-defined actions
584	CALL user_actions( 'before_timestep' )
585
586	!
587	!-- Calculate forces by wind turbines
588	IF ( wind_turbine ) THEN
589
590	CALL cpu_log( log_point(55), 'wind_turbine', 'start' )
591
592	CALL wtm_forces
593
594	CALL cpu_log( log_point(55), 'wind_turbine', 'stop' )
595
596	ENDIF
597
598	!
599	!-- Start of intermediate step loop
600	intermediate_timestep_count = 0
601	DO WHILE ( intermediate_timestep_count < &
602	intermediate_timestep_count_max )
603
604	intermediate_timestep_count = intermediate_timestep_count + 1
605
606	!
607	!-- Set the steering factors for the prognostic equations which depend
608	!-- on the timestep scheme
609	CALL timestep_scheme_steering
610
611	!
612	!-- Calculate those variables needed in the tendency terms which need
613	!-- global communication
614	IF ( .NOT. use_single_reference_value .AND. &
615	.NOT. use_initial_profile_as_reference ) THEN
616	!
617	!-- Horizontally averaged profiles to be used as reference state in
618	!-- buoyancy terms (WARNING: only the respective last call of
619	!-- calc_mean_profile defines the reference state!)
620	IF ( .NOT. neutral ) THEN
621	CALL calc_mean_profile( pt, 4 )
622	ref_state(:) = hom(:,1,4,0) ! this is used in the buoyancy term
623	ENDIF
624	IF ( ocean ) THEN
625	CALL calc_mean_profile( rho_ocean, 64 )
626	ref_state(:) = hom(:,1,64,0)
627	ENDIF
628	IF ( humidity ) THEN
629	CALL calc_mean_profile( vpt, 44 )
630	ref_state(:) = hom(:,1,44,0)
631	ENDIF
632	!
633	!-- Assure that ref_state does not become zero at any level
634	!-- ( might be the case if a vertical level is completely occupied
635	!-- with topography ).
636	ref_state = MERGE( MAXVAL(ref_state), ref_state, &
637	ref_state == 0.0_wp )
638	ENDIF
639
640	IF ( .NOT. constant_diffusion ) CALL production_e_init
641	IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. &
642	intermediate_timestep_count == 1 ) CALL ws_statistics
643	!
644	!-- In case of nudging calculate current nudging time scale and horizontal
645	!-- means of u, v, pt and q
646	IF ( nudging ) THEN
647	CALL calc_tnudge( simulated_time )
648	CALL calc_mean_profile( u, 1 )
649	CALL calc_mean_profile( v, 2 )
650	CALL calc_mean_profile( pt, 4 )
651	CALL calc_mean_profile( q, 41 )
652	ENDIF
653
654	!
655	!-- Solve the prognostic equations. A fast cache optimized version with
656	!-- only one single loop is used in case of Piascek-Williams advection
657	!-- scheme. NEC vector machines use a different version, because
658	!-- in the other versions a good vectorization is prohibited due to
659	!-- inlining problems.
660	IF ( loop_optimization == 'cache' ) THEN
661	CALL prognostic_equations_cache
662	ELSEIF ( loop_optimization == 'vector' ) THEN
663	CALL prognostic_equations_vector
664	ENDIF
665
666	!
667	!-- Particle transport/physics with the Lagrangian particle model
668	!-- (only once during intermediate steps, because it uses an Euler-step)
669	!-- ### particle model should be moved before prognostic_equations, in order
670	!-- to regard droplet interactions directly
671	IF ( particle_advection .AND. &
672	simulated_time >= particle_advection_start .AND. &
673	intermediate_timestep_count == 1 ) THEN
674	CALL lpm
675	first_call_lpm = .FALSE.
676	ENDIF
677
678	!
679	!-- Interaction of droplets with temperature and specific humidity.
680	!-- Droplet condensation and evaporation is calculated within
681	!-- advec_particles.
682	IF ( cloud_droplets .AND. &
683	intermediate_timestep_count == intermediate_timestep_count_max )&
684	THEN
685	CALL interaction_droplets_ptq
686	ENDIF
687
688	!
689	!-- Exchange of ghost points (lateral boundary conditions)
690	CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'start' )
691
692	CALL exchange_horiz( u_p, nbgp )
693	CALL exchange_horiz( v_p, nbgp )
694	CALL exchange_horiz( w_p, nbgp )
695	CALL exchange_horiz( pt_p, nbgp )
696	IF ( .NOT. constant_diffusion ) CALL exchange_horiz( e_p, nbgp )
697	IF ( rans_tke_e .OR. wang_kernel .OR. collision_turbulence &
698	.OR. use_sgs_for_particles ) THEN
699	IF ( rans_tke_e ) THEN
700	CALL exchange_horiz( diss_p, nbgp )
701	ELSE
702	CALL exchange_horiz( diss, nbgp )
703	ENDIF
704	ENDIF
705	IF ( ocean ) THEN
706	CALL exchange_horiz( sa_p, nbgp )
707	CALL exchange_horiz( rho_ocean, nbgp )
708	CALL exchange_horiz( prho, nbgp )
709	ENDIF
710	IF ( humidity ) THEN
711	CALL exchange_horiz( q_p, nbgp )
712	IF ( cloud_physics .AND. microphysics_morrison ) THEN
713	CALL exchange_horiz( qc_p, nbgp )
714	CALL exchange_horiz( nc_p, nbgp )
715	ENDIF
716	IF ( cloud_physics .AND. microphysics_seifert ) THEN
717	CALL exchange_horiz( qr_p, nbgp )
718	CALL exchange_horiz( nr_p, nbgp )
719	ENDIF
720	ENDIF
721	IF ( cloud_droplets ) THEN
722	CALL exchange_horiz( ql, nbgp )
723	CALL exchange_horiz( ql_c, nbgp )
724	CALL exchange_horiz( ql_v, nbgp )
725	CALL exchange_horiz( ql_vp, nbgp )
726	ENDIF
727	IF ( passive_scalar ) CALL exchange_horiz( s_p, nbgp )
728	IF ( air_chemistry ) THEN
729	DO n = 1, nspec
730	CALL exchange_horiz( chem_species(n)%conc_p, nbgp )
731	ENDDO
732	ENDIF
733
734	CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
735
736	!
737	!-- Boundary conditions for the prognostic quantities (except of the
738	!-- velocities at the outflow in case of a non-cyclic lateral wall)
739	CALL boundary_conds
740	!
741	!-- Swap the time levels in preparation for the next time step.
742	CALL swap_timelevel
743
744	!
745	!-- Vertical nesting: Interpolate fine grid data to the coarse grid
746	IF ( vnest_init ) THEN
747	CALL cpu_log( log_point(81), 'vnest_anterpolate', 'start' )
748	CALL vnest_anterpolate
749	CALL cpu_log( log_point(81), 'vnest_anterpolate', 'stop' )
750	ENDIF
751
752	IF ( nested_run ) THEN
753
754	CALL cpu_log( log_point(60), 'nesting', 'start' )
755	!
756	!-- Domain nesting. The data transfer subroutines pmci_parent_datatrans
757	!-- and pmci_child_datatrans are called inside the wrapper
758	!-- subroutine pmci_datatrans according to the control parameters
759	!-- nesting_mode and nesting_datatransfer_mode.
760	!-- TO_DO: why is nesting_mode given as a parameter here?
761	CALL pmci_datatrans( nesting_mode )
762
763	IF ( TRIM( nesting_mode ) == 'two-way' .OR. &
764	nesting_mode == 'vertical' ) THEN
765	!
766	!-- Exchange_horiz is needed for all parent-domains after the
767	!-- anterpolation
768	CALL exchange_horiz( u, nbgp )
769	CALL exchange_horiz( v, nbgp )
770	CALL exchange_horiz( w, nbgp )
771	IF ( .NOT. neutral ) CALL exchange_horiz( pt, nbgp )
772
773	IF ( humidity ) THEN
774
775	CALL exchange_horiz( q, nbgp )
776
777	IF ( cloud_physics .AND. microphysics_morrison ) THEN
778	CALL exchange_horiz( qc, nbgp )
779	CALL exchange_horiz( nc, nbgp )
780	ENDIF
781	IF ( cloud_physics .AND. microphysics_seifert ) THEN
782	CALL exchange_horiz( qr, nbgp )
783	CALL exchange_horiz( nr, nbgp )
784	ENDIF
785
786	ENDIF
787
788	IF ( passive_scalar ) CALL exchange_horiz( s, nbgp )
789	IF ( .NOT. constant_diffusion ) CALL exchange_horiz( e, nbgp )
790
791	IF ( .NOT. constant_diffusion .AND. rans_mode .AND. &
792	rans_tke_e ) &
793	CALL exchange_horiz( diss, nbgp )
794
795	IF ( air_chemistry ) THEN
796	DO n = 1, nspec
797	CALL exchange_horiz( chem_species(n)%conc, nbgp )
798	ENDDO
799	ENDIF
800
801	ENDIF
802	!
803	!-- Set boundary conditions again after interpolation and anterpolation.
804	CALL pmci_boundary_conds
805	!
806	!-- Correct the w top-BC in nest domains to ensure mass conservation.
807	!-- This action must never be done for the root domain. Vertical
808	!-- nesting implies mass conservation.
809	IF ( nest_domain ) THEN
810	CALL pmci_ensure_nest_mass_conservation
811	ENDIF
812
813	CALL cpu_log( log_point(60), 'nesting', 'stop' )
814
815	ENDIF
816
817	!
818	!-- Temperature offset must be imposed at cyclic boundaries in x-direction
819	!-- when a sloping surface is used
820	IF ( sloping_surface ) THEN
821	IF ( nxl == 0 ) pt(:,:,nxlg:nxl-1) = pt(:,:,nxlg:nxl-1) - &
822	pt_slope_offset
823	IF ( nxr == nx ) pt(:,:,nxr+1:nxrg) = pt(:,:,nxr+1:nxrg) + &
824	pt_slope_offset
825	ENDIF
826
827	!
828	!-- Impose a turbulent inflow using the recycling method
829	IF ( turbulent_inflow ) CALL inflow_turbulence
830
831	!
832	!-- Impose a turbulent inflow using synthetic generated turbulence
833	IF ( use_syn_turb_gen ) CALL stg_main
834
835	!
836	!-- Set values at outflow boundary using the special outflow condition
837	IF ( turbulent_outflow ) CALL outflow_turbulence
838
839	!
840	!-- Impose a random perturbation on the horizontal velocity field
841	IF ( create_disturbances .AND. &
842	( call_psolver_at_all_substeps .AND. &
843	intermediate_timestep_count == intermediate_timestep_count_max )&
844	.OR. ( .NOT. call_psolver_at_all_substeps .AND. &
845	intermediate_timestep_count == 1 ) ) &
846	THEN
847	time_disturb = time_disturb + dt_3d
848	IF ( time_disturb >= dt_disturb ) THEN
849	IF ( disturbance_energy_limit /= 0.0_wp .AND. &
850	hom(nzb+5,1,pr_palm,0) < disturbance_energy_limit ) THEN
851	CALL disturb_field( 'u', tend, u )
852	CALL disturb_field( 'v', tend, v )
853	ELSEIF ( ( .NOT. bc_lr_cyc .OR. .NOT. bc_ns_cyc ) &
854	.AND. .NOT. nest_domain .AND. .NOT. forcing ) THEN
855	!
856	!-- Runs with a non-cyclic lateral wall need perturbations
857	!-- near the inflow throughout the whole simulation
858	dist_range = 1
859	CALL disturb_field( 'u', tend, u )
860	CALL disturb_field( 'v', tend, v )
861	dist_range = 0
862	ENDIF
863	time_disturb = time_disturb - dt_disturb
864	ENDIF
865	ENDIF
866
867	!
868	!-- Map forcing data derived from larger scale model onto domain
869	!-- boundaries.
870	IF ( forcing .AND. intermediate_timestep_count == &
871	intermediate_timestep_count_max ) THEN
872	CALL forcing_bc
873	!
874	!-- Moreover, ensure mass conservation
875	CALL forcing_bc_mass_conservation
876	ENDIF
877
878	!
879	!-- Reduce the velocity divergence via the equation for perturbation
880	!-- pressure.
881	IF ( intermediate_timestep_count == 1 .OR. &
882	call_psolver_at_all_substeps ) THEN
883
884	IF ( vnest_init ) THEN
885	!
886	!-- Compute pressure in the CG, interpolate top boundary conditions
887	!-- to the FG and then compute pressure in the FG
888	IF ( coupling_mode == 'vnested_crse' ) CALL pres
889
890	CALL cpu_log( log_point(82), 'vnest_bc', 'start' )
891	CALL vnest_boundary_conds
892	CALL cpu_log( log_point(82), 'vnest_bc', 'stop' )
893
894	IF ( coupling_mode == 'vnested_fine' ) CALL pres
895
896	!-- Anterpolate TKE, satisfy Germano Identity
897	CALL cpu_log( log_point(83), 'vnest_anter_e', 'start' )
898	CALL vnest_anterpolate_e
899	CALL cpu_log( log_point(83), 'vnest_anter_e', 'stop' )
900
901	ELSE
902
903	CALL pres
904
905	ENDIF
906
907	ENDIF
908
909	!
910	!-- If required, compute liquid water content
911	IF ( cloud_physics ) THEN
912	CALL calc_liquid_water_content
913	ENDIF
914	!
915	!-- If required, compute virtual potential temperature
916	IF ( humidity ) THEN
917	CALL compute_vpt
918	ENDIF
919
920	!
921	!-- Compute the diffusion quantities
922	IF ( .NOT. constant_diffusion ) THEN
923
924	!
925	!-- Determine surface fluxes shf and qsws and surface values
926	!-- pt_surface and q_surface in dependence on data from external
927	!-- file LSF_DATA respectively
928	IF ( ( large_scale_forcing .AND. lsf_surf ) .AND. &
929	intermediate_timestep_count == intermediate_timestep_count_max )&
930	THEN
931	CALL ls_forcing_surf( simulated_time )
932	ENDIF
933
934	!
935	!-- First the vertical (and horizontal) fluxes in the surface
936	!-- (constant flux) layer are computed
937	IF ( constant_flux_layer ) THEN
938	CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'start' )
939	CALL surface_layer_fluxes
940	CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'stop' )
941	ENDIF
942	!
943	!-- If required, solve the energy balance for the surface and run soil
944	!-- model. Call for horizontal as well as vertical surfaces
945	IF ( land_surface .AND. time_since_reference_point >= skip_time_do_lsm) THEN
946
947	CALL cpu_log( log_point(54), 'land_surface', 'start' )
948	!
949	!-- Call for horizontal upward-facing surfaces
950	CALL lsm_energy_balance( .TRUE., -1 )
951	CALL lsm_soil_model( .TRUE., -1, .TRUE. )
952	!
953	!-- Call for northward-facing surfaces
954	CALL lsm_energy_balance( .FALSE., 0 )
955	CALL lsm_soil_model( .FALSE., 0, .TRUE. )
956	!
957	!-- Call for southward-facing surfaces
958	CALL lsm_energy_balance( .FALSE., 1 )
959	CALL lsm_soil_model( .FALSE., 1, .TRUE. )
960	!
961	!-- Call for eastward-facing surfaces
962	CALL lsm_energy_balance( .FALSE., 2 )
963	CALL lsm_soil_model( .FALSE., 2, .TRUE. )
964	!
965	!-- Call for westward-facing surfaces
966	CALL lsm_energy_balance( .FALSE., 3 )
967	CALL lsm_soil_model( .FALSE., 3, .TRUE. )
968	!
969	!-- At the end, set boundary conditons for potential temperature
970	!-- and humidity after running the land-surface model. This
971	!-- might be important for the nesting, where arrays are transfered.
972	CALL lsm_boundary_condition
973
974	CALL cpu_log( log_point(54), 'land_surface', 'stop' )
975	ENDIF
976	!
977	!-- If required, solve the energy balance for urban surfaces and run
978	!-- the material heat model
979	IF (urban_surface) THEN
980	CALL cpu_log( log_point(74), 'urban_surface', 'start' )
981	CALL usm_surface_energy_balance
982	IF ( usm_material_model ) THEN
983	CALL usm_green_heat_model
984	CALL usm_material_heat_model
985	ENDIF
986
987	CALL usm_temperature_near_surface
988	!
989	!-- At the end, set boundary conditons for potential temperature
990	!-- and humidity after running the urban-surface model. This
991	!-- might be important for the nesting, where arrays are transfered.
992	CALL usm_boundary_condition
993
994	CALL cpu_log( log_point(74), 'urban_surface', 'stop' )
995	ENDIF
996	!
997	!-- Compute the diffusion coefficients
998	CALL cpu_log( log_point(17), 'diffusivities', 'start' )
999	IF ( .NOT. humidity ) THEN
1000	IF ( ocean ) THEN
1001	CALL tcm_diffusivities( prho, prho_reference )
1002	ELSE
1003	CALL tcm_diffusivities( pt, pt_reference )
1004	ENDIF
1005	ELSE
1006	CALL tcm_diffusivities( vpt, pt_reference )
1007	ENDIF
1008	CALL cpu_log( log_point(17), 'diffusivities', 'stop' )
1009	!
1010	!-- Vertical nesting: set fine grid eddy viscosity top boundary condition
1011	IF ( vnest_init ) CALL vnest_boundary_conds_khkm
1012
1013	ENDIF
1014
1015	!
1016	!-- If required, calculate radiative fluxes and heating rates
1017	IF ( radiation .AND. intermediate_timestep_count &
1018	== intermediate_timestep_count_max .AND. time_since_reference_point > &
1019	skip_time_do_radiation ) THEN
1020
1021	time_radiation = time_radiation + dt_3d
1022
1023	IF ( time_radiation >= dt_radiation .OR. force_radiation_call ) &
1024	THEN
1025
1026	CALL cpu_log( log_point(50), 'radiation', 'start' )
1027
1028	IF ( .NOT. force_radiation_call ) THEN
1029	time_radiation = time_radiation - dt_radiation
1030	ENDIF
1031
1032	CALL radiation_control
1033
1034	CALL cpu_log( log_point(50), 'radiation', 'stop' )
1035
1036	IF ( urban_surface .OR. land_surface .AND. &
1037	radiation_interactions ) THEN
1038	CALL cpu_log( log_point(75), 'radiation_interaction', 'start' )
1039	CALL radiation_interaction
1040	CALL cpu_log( log_point(75), 'radiation_interaction', 'stop' )
1041	ENDIF
1042
1043	ENDIF
1044	ENDIF
1045
1046	ENDDO ! Intermediate step loop
1047	!
1048	!-- If required, consider chemical emissions
1049	!-- (todo (FK): Implement hourly call of emissions, using time_utc from
1050	!-- data_and_time_mod.f90;
1051	!-- move the CALL to appropriate location)
1052	IF ( air_chemistry ) THEN
1053	CALL chem_emissions
1054	ENDIF
1055	!
1056	!-- If required, do UV exposure calculations
1057	IF ( uv_exposure ) THEN
1058	CALL uvem_calc_exposure
1059	ENDIF
1060	!
1061	!-- Increase simulation time and output times
1062	nr_timesteps_this_run = nr_timesteps_this_run + 1
1063	current_timestep_number = current_timestep_number + 1
1064	simulated_time = simulated_time + dt_3d
1065	time_since_reference_point = simulated_time - coupling_start_time
1066	simulated_time_chr = time_to_string( time_since_reference_point )
1067
1068
1069
1070
1071	IF ( simulated_time >= skip_time_data_output_av ) THEN
1072	time_do_av = time_do_av + dt_3d
1073	ENDIF
1074	IF ( simulated_time >= skip_time_do2d_xy ) THEN
1075	time_do2d_xy = time_do2d_xy + dt_3d
1076	ENDIF
1077	IF ( simulated_time >= skip_time_do2d_xz ) THEN
1078	time_do2d_xz = time_do2d_xz + dt_3d
1079	ENDIF
1080	IF ( simulated_time >= skip_time_do2d_yz ) THEN
1081	time_do2d_yz = time_do2d_yz + dt_3d
1082	ENDIF
1083	IF ( simulated_time >= skip_time_do3d ) THEN
1084	time_do3d = time_do3d + dt_3d
1085	ENDIF
1086	DO mid = 1, masks
1087	IF ( simulated_time >= skip_time_domask(mid) ) THEN
1088	time_domask(mid)= time_domask(mid) + dt_3d
1089	ENDIF
1090	ENDDO
1091	time_dvrp = time_dvrp + dt_3d
1092	IF ( simulated_time >= skip_time_dosp ) THEN
1093	time_dosp = time_dosp + dt_3d
1094	ENDIF
1095	time_dots = time_dots + dt_3d
1096	IF ( .NOT. first_call_lpm ) THEN
1097	time_dopts = time_dopts + dt_3d
1098	ENDIF
1099	IF ( simulated_time >= skip_time_dopr ) THEN
1100	time_dopr = time_dopr + dt_3d
1101	ENDIF
1102	time_dopr_listing = time_dopr_listing + dt_3d
1103	time_run_control = time_run_control + dt_3d
1104
1105	!
1106	!-- Data exchange between coupled models
1107	IF ( coupling_mode /= 'uncoupled' .AND. run_coupled &
1108	.AND. .NOT. vnested ) THEN
1109	time_coupling = time_coupling + dt_3d
1110
1111	!
1112	!-- In case of model termination initiated by the local model
1113	!-- (terminate_coupled > 0), the coupler must be skipped because it would
1114	!-- cause an MPI intercomminucation hang.
1115	!-- If necessary, the coupler will be called at the beginning of the
1116	!-- next restart run.
1117	DO WHILE ( time_coupling >= dt_coupling .AND. terminate_coupled == 0 )
1118	CALL surface_coupler
1119	time_coupling = time_coupling - dt_coupling
1120	ENDDO
1121	ENDIF
1122
1123	!
1124	!-- Execute the gust module actions
1125	IF ( gust_module_enabled ) THEN
1126	CALL gust_actions( 'after_integration' )
1127	ENDIF
1128
1129	!
1130	!-- Execute user-defined actions
1131	CALL user_actions( 'after_integration' )
1132
1133	!
1134	!-- If Galilei transformation is used, determine the distance that the
1135	!-- model has moved so far
1136	IF ( galilei_transformation ) THEN
1137	advected_distance_x = advected_distance_x + u_gtrans * dt_3d
1138	advected_distance_y = advected_distance_y + v_gtrans * dt_3d
1139	ENDIF
1140
1141	!
1142	!-- Check, if restart is necessary (because cpu-time is expiring or
1143	!-- because it is forced by user) and set stop flag
1144	!-- This call is skipped if the remote model has already initiated a restart.
1145	IF ( .NOT. terminate_run ) CALL check_for_restart
1146
1147	!
1148	!-- Carry out statistical analysis and output at the requested output times.
1149	!-- The MOD function is used for calculating the output time counters (like
1150	!-- time_dopr) in order to regard a possible decrease of the output time
1151	!-- interval in case of restart runs
1152
1153	!
1154	!-- Set a flag indicating that so far no statistics have been created
1155	!-- for this time step
1156	flow_statistics_called = .FALSE.
1157
1158	!
1159	!-- If required, call flow_statistics for averaging in time
1160	IF ( averaging_interval_pr /= 0.0_wp .AND. &
1161	( dt_dopr - time_dopr ) <= averaging_interval_pr .AND. &
1162	simulated_time >= skip_time_dopr ) THEN
1163	time_dopr_av = time_dopr_av + dt_3d
1164	IF ( time_dopr_av >= dt_averaging_input_pr ) THEN
1165	do_sum = .TRUE.
1166	time_dopr_av = MOD( time_dopr_av, &
1167	MAX( dt_averaging_input_pr, dt_3d ) )
1168	ENDIF
1169	ENDIF
1170	IF ( do_sum ) CALL flow_statistics
1171
1172	!
1173	!-- Sum-up 3d-arrays for later output of time-averaged 2d/3d/masked data
1174	IF ( averaging_interval /= 0.0_wp .AND. &
1175	( dt_data_output_av - time_do_av ) <= averaging_interval .AND. &
1176	simulated_time >= skip_time_data_output_av ) &
1177	THEN
1178	time_do_sla = time_do_sla + dt_3d
1179	IF ( time_do_sla >= dt_averaging_input ) THEN
1180	CALL sum_up_3d_data
1181	average_count_3d = average_count_3d + 1
1182	time_do_sla = MOD( time_do_sla, MAX( dt_averaging_input, dt_3d ) )
1183	ENDIF
1184	ENDIF
1185
1186	!
1187	!-- Calculate spectra for time averaging
1188	IF ( averaging_interval_sp /= 0.0_wp .AND. &
1189	( dt_dosp - time_dosp ) <= averaging_interval_sp .AND. &
1190	simulated_time >= skip_time_dosp ) THEN
1191	time_dosp_av = time_dosp_av + dt_3d
1192	IF ( time_dosp_av >= dt_averaging_input_pr ) THEN
1193	CALL calc_spectra
1194	time_dosp_av = MOD( time_dosp_av, &
1195	MAX( dt_averaging_input_pr, dt_3d ) )
1196	ENDIF
1197	ENDIF
1198
1199	!
1200	!-- Call flight module and output data
1201	IF ( virtual_flight ) THEN
1202	CALL flight_measurement
1203	CALL data_output_flight
1204	ENDIF
1205
1206	!
1207	!-- Profile output (ASCII) on file
1208	IF ( time_dopr_listing >= dt_dopr_listing ) THEN
1209	CALL print_1d
1210	time_dopr_listing = MOD( time_dopr_listing, MAX( dt_dopr_listing, &
1211	dt_3d ) )
1212	ENDIF
1213
1214	!
1215	!-- Graphic output for PROFIL
1216	IF ( time_dopr >= dt_dopr ) THEN
1217	IF ( dopr_n /= 0 ) CALL data_output_profiles
1218	time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
1219	time_dopr_av = 0.0_wp ! due to averaging (see above)
1220	ENDIF
1221
1222	!
1223	!-- Graphic output for time series
1224	IF ( time_dots >= dt_dots ) THEN
1225	CALL data_output_tseries
1226	time_dots = MOD( time_dots, MAX( dt_dots, dt_3d ) )
1227	ENDIF
1228
1229	!
1230	!-- Output of spectra (formatted for use with PROFIL), in case of no
1231	!-- time averaging, spectra has to be calculated before
1232	IF ( time_dosp >= dt_dosp ) THEN
1233	IF ( average_count_sp == 0 ) CALL calc_spectra
1234	CALL data_output_spectra
1235	time_dosp = MOD( time_dosp, MAX( dt_dosp, dt_3d ) )
1236	ENDIF
1237
1238	!
1239	!-- 2d-data output (cross-sections)
1240	IF ( time_do2d_xy >= dt_do2d_xy ) THEN
1241	CALL data_output_2d( 'xy', 0 )
1242	time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
1243	ENDIF
1244	IF ( time_do2d_xz >= dt_do2d_xz ) THEN
1245	CALL data_output_2d( 'xz', 0 )
1246	time_do2d_xz = MOD( time_do2d_xz, MAX( dt_do2d_xz, dt_3d ) )
1247	ENDIF
1248	IF ( time_do2d_yz >= dt_do2d_yz ) THEN
1249	CALL data_output_2d( 'yz', 0 )
1250	time_do2d_yz = MOD( time_do2d_yz, MAX( dt_do2d_yz, dt_3d ) )
1251	ENDIF
1252
1253	!
1254	!-- 3d-data output (volume data)
1255	IF ( time_do3d >= dt_do3d ) THEN
1256	CALL data_output_3d( 0 )
1257	time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
1258	ENDIF
1259
1260	!
1261	!-- Masked data output
1262	DO mid = 1, masks
1263	IF ( time_domask(mid) >= dt_domask(mid) ) THEN
1264	CALL data_output_mask( 0 )
1265	time_domask(mid) = MOD( time_domask(mid), &
1266	MAX( dt_domask(mid), dt_3d ) )
1267	ENDIF
1268	ENDDO
1269
1270	!
1271	!-- Output of time-averaged 2d/3d/masked data
1272	IF ( time_do_av >= dt_data_output_av ) THEN
1273	CALL average_3d_data
1274	CALL data_output_2d( 'xy', 1 )
1275	CALL data_output_2d( 'xz', 1 )
1276	CALL data_output_2d( 'yz', 1 )
1277	CALL data_output_3d( 1 )
1278	DO mid = 1, masks
1279	CALL data_output_mask( 1 )
1280	ENDDO
1281	time_do_av = MOD( time_do_av, MAX( dt_data_output_av, dt_3d ) )
1282	ENDIF
1283
1284	!
1285	!-- Output of particle time series
1286	IF ( particle_advection ) THEN
1287	IF ( time_dopts >= dt_dopts .OR. &
1288	( simulated_time >= particle_advection_start .AND. &
1289	first_call_lpm ) ) THEN
1290	CALL data_output_ptseries
1291	time_dopts = MOD( time_dopts, MAX( dt_dopts, dt_3d ) )
1292	ENDIF
1293	ENDIF
1294
1295	!
1296	!-- Output of dvrp-graphics (isosurface, particles, slicer)
1297	#if defined( __dvrp_graphics )
1298	CALL DVRP_LOG_EVENT( -2, current_timestep_number-1 )
1299	#endif
1300	IF ( time_dvrp >= dt_dvrp ) THEN
1301	CALL data_output_dvrp
1302	time_dvrp = MOD( time_dvrp, MAX( dt_dvrp, dt_3d ) )
1303	ENDIF
1304	#if defined( __dvrp_graphics )
1305	CALL DVRP_LOG_EVENT( 2, current_timestep_number )
1306	#endif
1307
1308	!
1309	!-- If required, set the heat flux for the next time step at a random value
1310	IF ( constant_heatflux .AND. random_heatflux ) THEN
1311	IF ( surf_def_h(0)%ns >= 1 ) CALL disturb_heatflux( surf_def_h(0) )
1312	IF ( surf_lsm_h%ns >= 1 ) CALL disturb_heatflux( surf_lsm_h )
1313	IF ( surf_usm_h%ns >= 1 ) CALL disturb_heatflux( surf_usm_h )
1314	ENDIF
1315
1316	!
1317	!-- Execute the gust module actions
1318	IF ( gust_module_enabled ) THEN
1319	CALL gust_actions( 'after_timestep' )
1320	ENDIF
1321
1322	!
1323	!-- Execute user-defined actions
1324	CALL user_actions( 'after_timestep' )
1325
1326	!
1327	!-- Determine size of next time step. Save timestep dt_3d because it is
1328	!-- newly calculated in routine timestep, but required further below for
1329	!-- steering the run control output interval
1330	dt_3d_old = dt_3d
1331	CALL timestep
1332
1333	!
1334	!-- Synchronize the timestep in case of nested run.
1335	IF ( nested_run ) THEN
1336	!
1337	!-- Synchronize by unifying the time step.
1338	!-- Global minimum of all time-steps is used for all.
1339	CALL pmci_synchronize
1340	ENDIF
1341
1342	!
1343	!-- Computation and output of run control parameters.
1344	!-- This is also done whenever perturbations have been imposed
1345	IF ( time_run_control >= dt_run_control .OR. &
1346	timestep_scheme(1:5) /= 'runge' .OR. disturbance_created ) &
1347	THEN
1348	CALL run_control
1349	IF ( time_run_control >= dt_run_control ) THEN
1350	time_run_control = MOD( time_run_control, &
1351	MAX( dt_run_control, dt_3d_old ) )
1352	ENDIF
1353	ENDIF
1354
1355	!
1356	!-- Output elapsed simulated time in form of a progress bar on stdout
1357	IF ( myid == 0 ) CALL output_progress_bar
1358
1359	CALL cpu_log( log_point_s(10), 'timesteps', 'stop' )
1360
1361
1362	ENDDO ! time loop
1363
1364	!-- Vertical nesting: Deallocate variables initialized for vertical nesting
1365	IF ( vnest_init ) CALL vnest_deallocate
1366
1367	IF ( myid == 0 ) CALL finish_progress_bar
1368
1369	#if defined( __dvrp_graphics )
1370	CALL DVRP_LOG_EVENT( -2, current_timestep_number )
1371	#endif
1372
1373	CALL location_message( 'finished time-stepping', .TRUE. )
1374
1375	END SUBROUTINE time_integration

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |