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

Last change on this file since 2938 was 2938, checked in by suehring, 4 years ago

Nesting in RANS-LES and RANS-RANS mode enabled; synthetic turbulence generator at all lateral boundaries in nesting or non-cyclic forcing mode; revised Inifor initialization in nesting mode

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