source: palm/tags/release-5.0/SOURCE/time_integration.f90 @ 3984

Last change on this file since 3984 was 2696, checked in by kanani, 6 years ago

Merge of branch palm4u into trunk

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