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

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

Synchronize parent and child model after initialization and spinup phase; Check for consistent setting of spinup times in parent and child model; remove obsolete masking of tendency arrays during initialization

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