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

Last change on this file since 2776 was 2776, checked in by Giersch, 4 years ago

Skipping of module related restart data changed + adapting synthetic turbulence generator to current restart procedure

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