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

Last change on this file since 3183 was 3183, checked in by suehring, 6 years ago

last commit documented

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