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

Last change on this file since 2977 was 2977, checked in by kanani, 4 years ago

Fixes for radiative transfer model

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