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

Last change on this file since 3049 was 3042, checked in by schwenkel, 6 years ago

Changed the name specific humidity to mixing ratio

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