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

Last change on this file since 4047 was 4047, checked in by knoop, 2 years ago

Initial introduction of the dynamics module with only dynamics_swap_timelevel implemented

  • Property svn:keywords set to Id
File size: 80.3 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-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: time_integration.f90 4047 2019-06-21 18:58:09Z knoop $
27! Added remainings of swap_timelevel upon its dissolution
28!
29! 4043 2019-06-18 16:59:00Z schwenkel
30! Further LPM modularization
31!
32! 4039 2019-06-18 10:32:41Z suehring
33! Rename subroutines in module for diagnostic quantities
34!
35! 4029 2019-06-14 14:04:35Z raasch
36! exchange of ghost points and boundary conditions separated for chemical species and SALSA module,
37! bugfix: decycling of chemistry species after nesting data transfer
38!
39! 4022 2019-06-12 11:52:39Z suehring
40! Call synthetic turbulence generator at last RK3 substep right after boundary
41! conditions are updated in offline nesting in order to assure that
42! perturbations are always imposed
43!
44! 4017 2019-06-06 12:16:46Z schwenkel
45! Mass (volume) flux correction included to ensure global mass conservation for child domains.
46!
47! 3994 2019-05-22 18:08:09Z suehring
48! output of turbulence intensity added
49!
50! 3988 2019-05-22 11:32:37Z kanani
51! Implement steerable output interval for virtual measurements
52!
53! 3968 2019-05-13 11:04:01Z suehring
54! replace nspec_out with n_matched_vars
55!
56! 3929 2019-04-24 12:52:08Z banzhafs
57! Reverse changes back from revision 3878: use chem_boundary_conds instead of
58! chem_boundary_conds_decycle
59!
60!
61! 3885 2019-04-11 11:29:34Z kanani
62! Changes related to global restructuring of location messages and introduction
63! of additional debug messages
64!
65! 3879 2019-04-08 20:25:23Z knoop
66! Moved wtm_forces to module_interface_actions
67!
68! 3872 2019-04-08 15:03:06Z knoop
69! Modifications made for salsa:
70! - Call salsa_emission_update at each time step but do the checks within
71!   salsa_emission_update (i.e. skip_time_do_salsa >= time_since_reference_point
72!   and next_aero_emission_update <= time_since_reference_point ).
73! - Renamed nbins --> nbins_aerosol, ncc_tot --> ncomponents_mass and
74!   ngast --> ngases_salsa and loop indices b, c and sg to ib, ic and ig
75! - Apply nesting for salsa variables
76! - Removed cpu_log calls speciffic for salsa.
77!
78! 3833 2019-03-28 15:04:04Z forkel
79! added USE chem_gasphase_mod, replaced nspec by nspec since fixed compounds are not integrated
80!
81! 3820 2019-03-27 11:53:41Z forkel
82! renamed do_emiss to emissions_anthropogenic (ecc)
83!
84!
85! 3774 2019-03-04 10:52:49Z moh.hefny
86! rephrase if statement to avoid unallocated array in case of
87! nesting_offline is false (crashing during debug mode)
88!
89! 3761 2019-02-25 15:31:42Z raasch $
90! module section re-formatted and openacc required variables moved to separate section,
91! re-formatting to 100 char line width
92!
93! 3745 2019-02-15 18:57:56Z suehring
94! Call indoor model after first timestep
95!
96! 3744 2019-02-15 18:38:58Z suehring
97! - Moved call of bio_calculate_thermal_index_maps from biometeorology module to
98! time_integration to make sure averaged input is updated before calculating.
99!
100! 3739 2019-02-13 08:05:17Z dom_dwd_user
101! Removed everything related to "time_bio_results" as this is never used.
102!
103! 3724 2019-02-06 16:28:23Z kanani
104! Correct double-used log_point_s unit
105!
106! 3719 2019-02-06 13:10:18Z kanani
107! - removed wind_turbine cpu measurement, since same time is measured inside
108!   wtm_forces subroutine as special measures
109! - moved the numerous vnest cpulog to special measures
110! - extended radiation cpulog over entire radiation part,
111!   moved radiation_interactions cpulog to special measures
112! - moved some cpu_log calls to this routine for better overview
113!
114! 3705 2019-01-29 19:56:39Z suehring
115! Data output for virtual measurements added
116!
117! 3704 2019-01-29 19:51:41Z suehring
118! Rename subroutines for surface-data output
119!
120! 3647 2019-01-02 14:10:44Z kanani
121! Bugfix: add time_since_reference_point to IF clause for data_output calls
122! (otherwise skip_time_* values don't come into affect with dt_do* = 0.0).
123! Clean up indoor_model and biometeorology model call.
124!
125! 3646 2018-12-28 17:58:49Z kanani
126! Bugfix: use time_since_reference_point instead of simulated_time where
127! required (relevant when using wall/soil spinup)
128!
129! 3634 2018-12-18 12:31:28Z knoop
130! OpenACC port for SPEC
131!
132! 3597 2018-12-04 08:40:18Z maronga
133! Removed call to calculation of near air (10 cm) potential temperature (now in
134! surface layer fluxes)
135!
136! 3589 2018-11-30 15:09:51Z suehring
137! Move the control parameter "salsa" from salsa_mod to control_parameters
138! (M. Kurppa)
139!
140! 3582 2018-11-29 19:16:36Z suehring
141! dom_dwd_user, Schrempf:
142! Changes due to merge of uv exposure model into biometeorology_mod.
143!
144! 3525 2018-11-14 16:06:14Z kanani
145! Changes related to clean-up of biometeorology (dom_dwd_user)
146!
147! 3524 2018-11-14 13:36:44Z raasch
148! unused variables removed
149!
150! 3484 2018-11-02 14:41:25Z hellstea
151! pmci_ensure_nest_mass_conservation is premanently removed
152!
153! 3473 2018-10-30 20:50:15Z suehring
154! new module for virtual measurements introduced
155!
156! 3472 2018-10-30 20:43:50Z suehring
157! Add indoor model (kanani, srissman, tlang)
158!
159! 3467 2018-10-30 19:05:21Z suehring
160! Implementation of a new aerosol module salsa.
161!
162! 3448 2018-10-29 18:14:31Z kanani
163! Add biometeorology
164!
165! 3421 2018-10-24 18:39:32Z gronemeier
166! Surface data output
167!
168! 3418 2018-10-24 16:07:39Z kanani
169! call to material_heat_model now with check if spinup runs (rvtils)
170!
171! 3378 2018-10-19 12:34:59Z kanani
172! merge from radiation branch (r3362) into trunk
173! (moh.hefny):
174! Bugfix in the if statement to call radiation_interaction
175!
176! 3347 2018-10-15 14:21:08Z suehring
177! - offline nesting separated from large-scale forcing module
178! - changes for synthetic turbulence generator
179!
180! 3343 2018-10-15 10:38:52Z suehring
181! - Formatting, clean-up, comments (kanani)
182! - Added CALL to chem_emissions_setup (Russo)
183! - Code added for decycling chemistry (basit)
184!
185! 3294 2018-10-01 02:37:10Z raasch
186! changes concerning modularization of ocean option
187!
188! 3274 2018-09-24 15:42:55Z knoop
189! Modularization of all bulk cloud physics code components
190!
191! 3241 2018-09-12 15:02:00Z raasch
192! unused variables removed
193!
194! 3198 2018-08-15 09:23:10Z sward
195! Added multi_agent_system_end; defined start time for MAS relative to
196! time_since_reference_point
197!
198! 3183 2018-07-27 14:25:55Z suehring
199! Replace simulated_time by time_since_reference_point in COSMO nesting mode.
200! Rename subroutines and variables in COSMO nesting mode
201!
202! 3182 2018-07-27 13:36:03Z suehring
203! Added multi agent system
204!
205! 3042 2018-05-25 10:44:37Z schwenkel
206! Changed the name specific humidity to mixing ratio
207!
208! 3040 2018-05-25 10:22:08Z schwenkel
209! Fixed bug in IF statement
210! Ensure that the time when calling the radiation to be the time step of the
211! pre-calculated time when first calculate the positions of the sun
212!
213! 3004 2018-04-27 12:33:25Z Giersch
214! First call of flow_statistics has been removed. It is already called in
215! run_control itself
216!
217! 2984 2018-04-18 11:51:30Z hellstea
218! CALL pmci_ensure_nest_mass_conservation is removed (so far only commented out)
219! as seemingly unnecessary.
220!
221! 2941 2018-04-03 11:54:58Z kanani
222! Deduct spinup_time from RUN_CONTROL output of main 3d run
223! (use time_since_reference_point instead of simulated_time)
224!
225! 2938 2018-03-27 15:52:42Z suehring
226! Nesting of dissipation rate in case of RANS mode and TKE-e closure is applied
227!
228! 2936 2018-03-27 14:49:27Z suehring
229! Little formatting adjustment.
230!
231! 2817 2018-02-19 16:32:21Z knoop
232! Preliminary gust module interface implemented
233!
234! 2801 2018-02-14 16:01:55Z thiele
235! Changed lpm from subroutine to module.
236! Introduce particle transfer in nested models.
237!
238! 2776 2018-01-31 10:44:42Z Giersch
239! Variable use_synthetic_turbulence_generator has been abbreviated
240!
241! 2773 2018-01-30 14:12:54Z suehring
242! - Nesting for chemical species
243!
244! 2766 2018-01-22 17:17:47Z kanani
245! Removed preprocessor directive __chem
246!
247! 2718 2018-01-02 08:49:38Z maronga
248! Corrected "Former revisions" section
249!
250! 2696 2017-12-14 17:12:51Z kanani
251! - Change in file header (GPL part)
252! - Implementation of uv exposure model (FK)
253! - Moved vnest_boundary_conds_khkm from tcm_diffusivities to here (TG)
254! - renamed diffusivities to tcm_diffusivities (TG)
255! - implement prognostic equation for diss (TG)
256! - Moved/commented CALL to chem_emissions (FK)
257! - Added CALL to chem_emissions (FK)
258! - Implementation of chemistry module (FK)
259! - Calls for setting boundary conditions in USM and LSM (MS)
260! - Large-scale forcing with larger-scale models implemented (MS)
261! - Rename usm_radiation into radiation_interactions; merge with branch
262!   radiation (MS)
263! - added call for usm_green_heat_model for green building surfaces (RvT)
264! - added call for usm_temperature_near_surface for use in indoor model (RvT)
265!
266! 2617 2017-11-16 12:47:24Z suehring
267! Bugfix, assure that the reference state does not become zero.
268!
269! 2563 2017-10-19 15:36:10Z Giersch
270! Variable wind_turbine moved to module control_parameters
271!
272! 2365 2017-08-21 14:59:59Z kanani
273! Vertical grid nesting implemented (SadiqHuq)
274!
275! 2320 2017-07-21 12:47:43Z suehring
276! Set bottom boundary conditions after nesting interpolation and anterpolation
277!
278! 2299 2017-06-29 10:14:38Z maronga
279! Call of soil model adjusted
280!
281! 2292 2017-06-20 09:51:42Z schwenkel
282! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
283! includes two more prognostic equations for cloud drop concentration (nc) 
284! and cloud water content (qc).
285!
286! 2271 2017-06-09 12:34:55Z sward
287! Start timestep message changed
288!
289! 2259 2017-06-08 09:09:11Z gronemeier
290! Implemented synthetic turbulence generator
291!
292! 2233 2017-05-30 18:08:54Z suehring
293!
294! 2232 2017-05-30 17:47:52Z suehring
295! Adjustments to new topography and surface concept
296! Modify passed parameters for disturb_field
297!
298! 2178 2017-03-17 11:07:39Z hellstea
299! Setting perturbations at all times near inflow boundary is removed
300! in case of nested boundaries
301!
302! 2174 2017-03-13 08:18:57Z maronga
303! Added support for nesting with cloud microphysics
304!
305! 2118 2017-01-17 16:38:49Z raasch
306! OpenACC directives and related code removed
307!
308! 2050 2016-11-08 15:00:55Z gronemeier
309! Implement turbulent outflow condition
310!
311! 2031 2016-10-21 15:11:58Z knoop
312! renamed variable rho to rho_ocean
313!
314! 2011 2016-09-19 17:29:57Z kanani
315! Flag urban_surface is now defined in module control_parameters,
316! removed commented CALLs of global_min_max.
317!
318! 2007 2016-08-24 15:47:17Z kanani
319! Added CALLs for new urban surface model
320!
321! 2000 2016-08-20 18:09:15Z knoop
322! Forced header and separation lines into 80 columns
323!
324! 1976 2016-07-27 13:28:04Z maronga
325! Simplified calls to radiation model
326!
327! 1960 2016-07-12 16:34:24Z suehring
328! Separate humidity and passive scalar
329!
330! 1957 2016-07-07 10:43:48Z suehring
331! flight module added
332!
333! 1919 2016-05-27 14:51:23Z raasch
334! Initial version of purely vertical nesting introduced.
335!
336! 1918 2016-05-27 14:35:57Z raasch
337! determination of time step moved to the end of the time step loop,
338! the first time step is now always calculated before the time step loop (i.e.
339! also in case of restart runs)
340!
341! 1914 2016-05-26 14:44:07Z witha
342! Added call for wind turbine model
343!
344! 1878 2016-04-19 12:30:36Z hellstea
345! Synchronization for nested runs rewritten
346!
347! 1853 2016-04-11 09:00:35Z maronga
348! Adjusted for use with radiation_scheme = constant
349!
350! 1849 2016-04-08 11:33:18Z hoffmann
351! Adapted for modularization of microphysics
352!
353! 1833 2016-04-07 14:23:03Z raasch
354! spectrum renamed spectra_mod, spectra related variables moved to spectra_mod
355!
356! 1831 2016-04-07 13:15:51Z hoffmann
357! turbulence renamed collision_turbulence
358!
359! 1822 2016-04-07 07:49:42Z hoffmann
360! icloud_scheme replaced by microphysics_*
361!
362! 1808 2016-04-05 19:44:00Z raasch
363! output message in case unscheduled radiation calls removed
364!
365! 1797 2016-03-21 16:50:28Z raasch
366! introduction of different datatransfer modes
367!
368! 1791 2016-03-11 10:41:25Z raasch
369! call of pmci_update_new removed
370!
371! 1786 2016-03-08 05:49:27Z raasch
372! +module spectrum
373!
374! 1783 2016-03-06 18:36:17Z raasch
375! switch back of netcdf data format for mask output moved to the mask output
376! routine
377!
378! 1781 2016-03-03 15:12:23Z raasch
379! some pmc calls removed at the beginning (before timeloop),
380! pmc initialization moved to the main program
381!
382! 1764 2016-02-28 12:45:19Z raasch
383! PMC_ACTIVE flags removed,
384! bugfix: nest synchronization after first call of timestep
385!
386! 1762 2016-02-25 12:31:13Z hellstea
387! Introduction of nested domain feature
388!
389! 1736 2015-12-04 08:56:33Z raasch
390! no perturbations added to total domain if energy limit has been set zero
391!
392! 1691 2015-10-26 16:17:44Z maronga
393! Added option for spin-ups without land surface and radiation models. Moved calls
394! for radiation and lan surface schemes.
395!
396! 1682 2015-10-07 23:56:08Z knoop
397! Code annotations made doxygen readable
398!
399! 1671 2015-09-25 03:29:37Z raasch
400! bugfix: ghostpoint exchange for array diss in case that sgs velocities are used
401! for particles
402!
403! 1585 2015-04-30 07:05:52Z maronga
404! Moved call of radiation scheme. Added support for RRTM
405!
406! 1551 2015-03-03 14:18:16Z maronga
407! Added interface for different radiation schemes.
408!
409! 1496 2014-12-02 17:25:50Z maronga
410! Added calls for the land surface model and radiation scheme
411!
412! 1402 2014-05-09 14:25:13Z raasch
413! location messages modified
414!
415! 1384 2014-05-02 14:31:06Z raasch
416! location messages added
417!
418! 1380 2014-04-28 12:40:45Z heinze
419! CALL of nudge_ref added
420! bc_pt_t_val and bc_q_t_val are updated in case nudging is used
421!
422! 1365 2014-04-22 15:03:56Z boeske
423! Reset sums_ls_l to zero at each timestep
424! +sums_ls_l
425! Calculation of reference state (previously in subroutine calc_mean_profile)
426
427! 1342 2014-03-26 17:04:47Z kanani
428! REAL constants defined as wp-kind
429!
430! 1320 2014-03-20 08:40:49Z raasch
431! ONLY-attribute added to USE-statements,
432! kind-parameters added to all INTEGER and REAL declaration statements,
433! kinds are defined in new module kinds,
434! old module precision_kind is removed,
435! revision history before 2012 removed,
436! comment fields (!:) to be used for variable explanations added to
437! all variable declaration statements
438! 1318 2014-03-17 13:35:16Z raasch
439! module interfaces removed
440!
441! 1308 2014-03-13 14:58:42Z fricke
442! +netcdf_data_format_save
443! For masked data, parallel netcdf output is not tested so far, hence
444! netcdf_data_format is switched back to non-paralell output.
445!
446! 1276 2014-01-15 13:40:41Z heinze
447! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
448!
449! 1257 2013-11-08 15:18:40Z raasch
450! acc-update-host directive for timestep removed
451!
452! 1241 2013-10-30 11:36:58Z heinze
453! Generalize calc_mean_profile for wider use
454! Determine shf and qsws in dependence on data from LSF_DATA
455! Determine ug and vg in dependence on data from LSF_DATA
456! 1221 2013-09-10 08:59:13Z raasch
457! host update of arrays before timestep is called
458!
459! 1179 2013-06-14 05:57:58Z raasch
460! mean profiles for reference state are only calculated if required,
461! small bugfix for background communication
462!
463! 1171 2013-05-30 11:27:45Z raasch
464! split of prognostic_equations deactivated (comment lines), for the time being
465!
466! 1128 2013-04-12 06:19:32Z raasch
467! asynchronous transfer of ghost point data realized for acc-optimized version:
468! prognostic_equations are first called two times for those points required for
469! the left-right and north-south exchange, respectively, and then for the
470! remaining points,
471! those parts requiring global communication moved from prognostic_equations to
472! here
473!
474! 1115 2013-03-26 18:16:16Z hoffmann
475! calculation of qr and nr is restricted to precipitation
476!
477! 1113 2013-03-10 02:48:14Z raasch
478! GPU-porting of boundary conditions,
479! openACC directives updated
480! formal parameter removed from routine boundary_conds
481!
482! 1111 2013-03-08 23:54:10Z raasch
483! +internal timestep counter for cpu statistics added,
484! openACC directives updated
485!
486! 1092 2013-02-02 11:24:22Z raasch
487! unused variables removed
488!
489! 1065 2012-11-22 17:42:36Z hoffmann
490! exchange of diss (dissipation rate) in case of turbulence = .TRUE. added
491!
492! 1053 2012-11-13 17:11:03Z hoffmann
493! exchange of ghost points for nr, qr added
494!
495! 1036 2012-10-22 13:43:42Z raasch
496! code put under GPL (PALM 3.9)
497!
498! 1019 2012-09-28 06:46:45Z raasch
499! non-optimized version of prognostic_equations removed
500!
501! 1015 2012-09-27 09:23:24Z raasch
502! +call of prognostic_equations_acc
503!
504! 1001 2012-09-13 14:08:46Z raasch
505! all actions concerning leapfrog- and upstream-spline-scheme removed
506!
507! 849 2012-03-15 10:35:09Z raasch
508! advec_particles renamed lpm, first_call_advec_particles renamed first_call_lpm
509!
510! 825 2012-02-19 03:03:44Z raasch
511! wang_collision_kernel renamed wang_kernel
512!
513! Revision 1.1  1997/08/11 06:19:04  raasch
514! Initial revision
515!
516!
517! Description:
518! ------------
519!> Integration in time of the model equations, statistical analysis and graphic
520!> output
521!------------------------------------------------------------------------------!
522 SUBROUTINE time_integration
523 
524
525    USE advec_ws,                                                                                  &
526        ONLY:  ws_statistics
527
528    USE arrays_3d,                                                                                 &
529        ONLY:  diss, diss_p, dzu, e, e_p, nc, nc_p, nr, nr_p, prho, pt, pt_p, pt_init, q_init, q,  &
530               qc, qc_p, ql, ql_c, ql_v, ql_vp, qr, qr_p, q_p, ref_state, rho_ocean, s, s_p, sa_p, &
531               tend, u, u_p, v, vpt, v_p, w, w_p
532
533    USE biometeorology_mod,                                                                        &
534        ONLY:  bio_calculate_thermal_index_maps, thermal_comfort, uvem_calc_exposure, uv_exposure
535
536    USE bulk_cloud_model_mod,                                                                      &
537        ONLY: bulk_cloud_model, calc_liquid_water_content, collision_turbulence,                   &
538              microphysics_morrison, microphysics_seifert
539
540    USE calc_mean_profile_mod,                                                                     &
541        ONLY:  calc_mean_profile
542
543    USE chem_emissions_mod,                                                                        &
544        ONLY:  chem_emissions_setup
545
546    USE chem_gasphase_mod,                                                                         &
547        ONLY:  nvar
548
549    USE chem_modules,                                                                              &
550        ONLY:  bc_cs_t_val, chem_species, cs_name, emissions_anthropogenic, n_matched_vars
551
552    USE chemistry_model_mod,                                                                       &
553        ONLY:  chem_boundary_conds
554
555    USE control_parameters,                                                                        &
556        ONLY:  advected_distance_x, advected_distance_y, air_chemistry, average_count_3d,          &
557               averaging_interval, averaging_interval_pr, bc_lr_cyc, bc_ns_cyc, bc_pt_t_val,       &
558               bc_q_t_val, biometeorology, call_psolver_at_all_substeps,  child_domain,            &
559               cloud_droplets, constant_flux_layer, constant_heatflux, create_disturbances,        &
560               dopr_n, constant_diffusion, coupling_mode, coupling_start_time,                     &
561               current_timestep_number, disturbance_created, disturbance_energy_limit, dist_range, &
562               do_sum, dt_3d, dt_averaging_input, dt_averaging_input_pr, dt_coupling,              &
563               dt_data_output_av, dt_disturb, dt_do2d_xy, dt_do2d_xz, dt_do2d_yz, dt_do3d,         &
564               dt_domask,dt_dopts, dt_dopr, dt_dopr_listing, dt_dots, dt_run_control,              &
565               end_time, first_call_lpm, first_call_mas, galilei_transformation, humidity,         &
566               indoor_model, intermediate_timestep_count, intermediate_timestep_count_max,         &
567               land_surface, large_scale_forcing, loop_optimization, lsf_surf, lsf_vert, masks,    &
568               mid, multi_agent_system_end, multi_agent_system_start, nesting_offline, neutral,    &
569               nr_timesteps_this_run, nudging, ocean_mode, passive_scalar, pt_reference,           &
570               pt_slope_offset, random_heatflux, rans_mode, rans_tke_e, run_coupled, salsa,        &
571               simulated_time, simulated_time_chr, skip_time_do2d_xy, skip_time_do2d_xz,           &
572               skip_time_do2d_yz, skip_time_do3d, skip_time_domask, skip_time_dopr,                &
573               skip_time_data_output_av, sloping_surface, stop_dt, surface_output,                 &
574               terminate_coupled, terminate_run, timestep_scheme, time_coupling, time_do2d_xy,     &
575               time_do2d_xz, time_do2d_yz, time_do3d, time_domask, time_dopr, time_dopr_av,        &
576               time_dopr_listing, time_dopts, time_dosp, time_dosp_av, time_dots, time_do_av,      &
577               time_do_sla, time_disturb, time_run_control, time_since_reference_point,            &
578               turbulent_inflow, turbulent_outflow, urban_surface,                                 &
579               use_initial_profile_as_reference, use_single_reference_value, u_gtrans, v_gtrans,   &
580               virtual_flight, virtual_measurement, ws_scheme_mom, ws_scheme_sca, timestep_count
581
582    USE cpulog,                                                                                    &
583        ONLY:  cpu_log, log_point, log_point_s
584
585    USE date_and_time_mod,                                                                         &
586        ONLY:  calc_date_and_time, hour_call_emis, hour_of_year
587
588    USE diagnostic_output_quantities_mod,                                                          &
589        ONLY:  doq_calculate,                                                                      &
590               timestep_number_at_prev_calc
591
592    USE flight_mod,                                                                                &
593        ONLY:  flight_measurement
594
595    USE indices,                                                                                   &
596        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, nzb, nzt
597
598    USE indoor_model_mod,                                                                          &
599        ONLY:  dt_indoor, im_main_heatcool, time_indoor
600
601    USE interfaces
602
603    USE kinds
604
605    USE land_surface_model_mod,                                                                    &
606        ONLY:  lsm_boundary_condition, lsm_energy_balance, lsm_soil_model, skip_time_do_lsm
607
608    USE lagrangian_particle_model_mod,                                                             &
609        ONLY:  lpm_data_output_ptseries, lpm_interaction_droplets_ptq
610
611    USE lsf_nudging_mod,                                                                           &
612        ONLY:  calc_tnudge, ls_forcing_surf, ls_forcing_vert, nudge_ref
613
614    USE module_interface,                                                                          &
615        ONLY:  module_interface_actions, module_interface_swap_timelevel
616
617    USE multi_agent_system_mod,                                                                    &
618        ONLY:  agents_active, multi_agent_system
619
620    USE nesting_offl_mod,                                                                          &
621        ONLY:  nesting_offl_bc, nesting_offl_mass_conservation
622
623    USE netcdf_data_input_mod,                                                                     &
624        ONLY:  chem_emis, chem_emis_att, nest_offl, netcdf_data_input_offline_nesting
625
626    USE ocean_mod,                                                                                 &
627        ONLY:  prho_reference
628
629    USE particle_attributes,                                                                       &
630        ONLY:  particle_advection, particle_advection_start, use_sgs_for_particles, wang_kernel
631
632    USE pegrid
633
634    USE pmc_interface,                                                                             &
635        ONLY:  nested_run, nesting_mode, pmci_boundary_conds, pmci_datatrans, pmci_synchronize,    &
636        pmci_ensure_nest_mass_conservation, pmci_ensure_nest_mass_conservation_vertical,           &
637        pmci_set_swaplevel
638
639    USE progress_bar,                                                                              &
640        ONLY:  finish_progress_bar, output_progress_bar
641
642    USE prognostic_equations_mod,                                                                  &
643        ONLY:  prognostic_equations_cache, prognostic_equations_vector
644
645    USE radiation_model_mod,                                                                       &
646        ONLY: dt_radiation, force_radiation_call, radiation, radiation_control,                    &
647              radiation_interaction, radiation_interactions, skip_time_do_radiation, time_radiation
648
649    USE salsa_mod,                                                                                 &
650        ONLY: aerosol_number, aerosol_mass, bc_am_t_val, bc_an_t_val, bc_gt_t_val,                 &
651              nbins_aerosol, ncomponents_mass, ngases_salsa, salsa_boundary_conds,                 &
652              salsa_emission_update, salsa_gas, salsa_gases_from_chem, skip_time_do_salsa
653
654    USE spectra_mod,                                                                               &
655        ONLY: average_count_sp, averaging_interval_sp, calc_spectra, dt_dosp, skip_time_dosp
656
657    USE statistics,                                                                                &
658        ONLY:  flow_statistics_called, hom, pr_palm, sums_ls_l
659
660
661    USE surface_layer_fluxes_mod,                                                                  &
662        ONLY:  surface_layer_fluxes
663
664    USE surface_data_output_mod,                                                                   &
665        ONLY:  average_count_surf, averaging_interval_surf, dt_dosurf, dt_dosurf_av,               &
666               surface_data_output, surface_data_output_averaging, skip_time_dosurf,               &
667               skip_time_dosurf_av, time_dosurf, time_dosurf_av
668
669    USE surface_mod,                                                                               &
670        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
671
672    USE synthetic_turbulence_generator_mod,                                                        &
673        ONLY:  dt_stg_call, dt_stg_adjust, parametrize_inflow_turbulence, stg_adjust, stg_main,    &
674               time_stg_adjust, time_stg_call, use_syn_turb_gen
675
676    USE turbulence_closure_mod,                                                                    &
677        ONLY:  tcm_diffusivities, production_e_init
678
679    USE urban_surface_mod,                                                                         &
680        ONLY:  usm_boundary_condition, usm_material_heat_model, usm_material_model,                &
681               usm_surface_energy_balance, usm_green_heat_model
682
683    USE vertical_nesting_mod,                                                                      &
684        ONLY:  vnested, vnest_anterpolate, vnest_anterpolate_e, vnest_boundary_conds,              &
685               vnest_boundary_conds_khkm, vnest_deallocate, vnest_init, vnest_init_fine,           &
686               vnest_start_time
687
688    USE virtual_measurement_mod,                                                                   &
689        ONLY:  dt_virtual_measurement,                                                             &
690               time_virtual_measurement,                                                           &
691               vm_data_output,                                                                     &
692               vm_sampling,                                                                        &
693               vm_time_start
694
695
696#if defined( _OPENACC )
697    USE arrays_3d,                                                             &
698        ONLY:  d, dd2zu, ddzu, ddzw, drho_air, drho_air_zw, dzw, heatflux_output_conversion, kh,   &
699               km, momentumflux_output_conversion, p, ptdf_x, ptdf_y, rdf, rdf_sc, rho_air,        &
700               rho_air_zw, te_m, tpt_m, tu_m, tv_m, tw_m, ug, u_init, u_stokes_zu, vg, v_init,     &
701               v_stokes_zu, zu
702
703    USE control_parameters,                                                                        &
704        ONLY:  tsc
705
706    USE indices,                                                                                   &
707        ONLY:  advc_flags_1, advc_flags_2, nyn, nyng, nys, nysg, nz, nzb_max, wall_flags_0
708
709    USE statistics,                                                                                &
710        ONLY:  rmask, statistic_regions, sums_l, sums_l_l, sums_us2_ws_l,                          &
711               sums_wsus_ws_l, sums_vs2_ws_l, sums_wsvs_ws_l, sums_ws2_ws_l, sums_wspts_ws_l,      &
712               sums_wsqs_ws_l, sums_wssas_ws_l, sums_wsqcs_ws_l, sums_wsqrs_ws_l, sums_wsncs_ws_l, &
713               sums_wsnrs_ws_l, sums_wsss_ws_l, weight_substep, sums_salsa_ws_l
714
715    USE surface_mod,                                                                               &
716        ONLY:  bc_h, enter_surface_arrays, exit_surface_arrays
717#endif
718
719
720    IMPLICIT NONE
721
722    CHARACTER (LEN=9) ::  time_to_string   !<
723
724    INTEGER(iwp)      ::  ib        !< index for aerosol size bins
725    INTEGER(iwp)      ::  ic        !< index for aerosol mass bins
726    INTEGER(iwp)      ::  icc       !< additional index for aerosol mass bins
727    INTEGER(iwp)      ::  ig        !< index for salsa gases
728    INTEGER(iwp)      ::  lsp
729    INTEGER(iwp)      ::  lsp_usr   !<
730    INTEGER(iwp)      ::  n         !< loop counter for chemistry species
731
732    REAL(wp) ::  dt_3d_old  !< temporary storage of timestep to be used for
733                            !< steering of run control output interval
734    REAL(wp) ::  time_since_reference_point_save  !< original value of
735                                                  !< time_since_reference_point
736
737
738!
739!-- Copy data from arrays_3d
740!$ACC DATA &
741!$ACC COPY(d(nzb+1:nzt,nys:nyn,nxl:nxr)) &
742!$ACC COPY(e(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
743!$ACC COPY(u(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
744!$ACC COPY(v(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
745!$ACC COPY(w(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
746!$ACC COPY(kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
747!$ACC COPY(km(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
748!$ACC COPY(p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
749!$ACC COPY(pt(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
750
751!$ACC DATA &
752!$ACC COPY(e_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
753!$ACC COPY(u_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
754!$ACC COPY(v_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
755!$ACC COPY(w_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
756!$ACC COPY(pt_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
757!$ACC COPY(tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
758!$ACC COPY(te_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
759!$ACC COPY(tu_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
760!$ACC COPY(tv_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
761!$ACC COPY(tw_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
762!$ACC COPY(tpt_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
763
764!$ACC DATA &
765!$ACC COPYIN(rho_air(nzb:nzt+1), drho_air(nzb:nzt+1)) &
766!$ACC COPYIN(rho_air_zw(nzb:nzt+1), drho_air_zw(nzb:nzt+1)) &
767!$ACC COPYIN(zu(nzb:nzt+1)) &
768!$ACC COPYIN(dzu(1:nzt+1), dzw(1:nzt+1)) &
769!$ACC COPYIN(ddzu(1:nzt+1), dd2zu(1:nzt)) &
770!$ACC COPYIN(ddzw(1:nzt+1)) &
771!$ACC COPYIN(heatflux_output_conversion(nzb:nzt+1)) &
772!$ACC COPYIN(momentumflux_output_conversion(nzb:nzt+1)) &
773!$ACC COPYIN(rdf(nzb+1:nzt), rdf_sc(nzb+1:nzt)) &
774!$ACC COPYIN(ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng)) &
775!$ACC COPYIN(ref_state(0:nz+1)) &
776!$ACC COPYIN(u_init(0:nz+1), v_init(0:nz+1)) &
777!$ACC COPYIN(u_stokes_zu(nzb:nzt+1), v_stokes_zu(nzb:nzt+1)) &
778!$ACC COPYIN(pt_init(0:nz+1)) &
779!$ACC COPYIN(ug(0:nz+1), vg(0:nz+1))
780
781!
782!-- Copy data from control_parameters
783!$ACC DATA &
784!$ACC COPYIN(tsc(1:5))
785
786!
787!-- Copy data from indices
788!$ACC DATA &
789!$ACC COPYIN(advc_flags_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
790!$ACC COPYIN(advc_flags_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
791!$ACC COPYIN(wall_flags_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
792
793!
794!-- Copy data from surface_mod
795!$ACC DATA &
796!$ACC COPYIN(bc_h(0:1)) &
797!$ACC COPYIN(bc_h(0)%i(1:bc_h(0)%ns)) &
798!$ACC COPYIN(bc_h(0)%j(1:bc_h(0)%ns)) &
799!$ACC COPYIN(bc_h(0)%k(1:bc_h(0)%ns)) &
800!$ACC COPYIN(bc_h(1)%i(1:bc_h(1)%ns)) &
801!$ACC COPYIN(bc_h(1)%j(1:bc_h(1)%ns)) &
802!$ACC COPYIN(bc_h(1)%k(1:bc_h(1)%ns))
803
804!
805!-- Copy data from statistics
806!$ACC DATA &
807!$ACC COPYIN(hom(0:nz+1,1:2,1:4,0)) &
808!$ACC COPYIN(rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions)) &
809!$ACC COPYIN(weight_substep(1:intermediate_timestep_count_max)) &
810!$ACC COPY(sums_l(nzb:nzt+1,1:pr_palm,0)) &
811!$ACC COPY(sums_l_l(nzb:nzt+1,0:statistic_regions,0)) &
812!$ACC COPY(sums_us2_ws_l(nzb:nzt+1,0)) &
813!$ACC COPY(sums_wsus_ws_l(nzb:nzt+1,0)) &
814!$ACC COPY(sums_vs2_ws_l(nzb:nzt+1,0)) &
815!$ACC COPY(sums_wsvs_ws_l(nzb:nzt+1,0)) &
816!$ACC COPY(sums_ws2_ws_l(nzb:nzt+1,0)) &
817!$ACC COPY(sums_wspts_ws_l(nzb:nzt+1,0)) &
818!$ACC COPY(sums_wssas_ws_l(nzb:nzt+1,0)) &
819!$ACC COPY(sums_wsqs_ws_l(nzb:nzt+1,0)) &
820!$ACC COPY(sums_wsqcs_ws_l(nzb:nzt+1,0)) &
821!$ACC COPY(sums_wsqrs_ws_l(nzb:nzt+1,0)) &
822!$ACC COPY(sums_wsncs_ws_l(nzb:nzt+1,0)) &
823!$ACC COPY(sums_wsnrs_ws_l(nzb:nzt+1,0)) &
824!$ACC COPY(sums_wsss_ws_l(nzb:nzt+1,0)) &
825!$ACC COPY(sums_salsa_ws_l(nzb:nzt+1,0))
826
827#if defined( _OPENACC )
828    CALL enter_surface_arrays
829#endif
830
831!
832!-- At beginning determine the first time step
833    CALL timestep
834!
835!-- Synchronize the timestep in case of nested run.
836    IF ( nested_run )  THEN
837!
838!--    Synchronization by unifying the time step.
839!--    Global minimum of all time-steps is used for all.
840       CALL pmci_synchronize
841    ENDIF
842
843!
844!-- Determine and print out the run control quantities before the first time
845!-- step of this run. For the initial run, some statistics (e.g. divergence)
846!-- need to be determined first --> CALL flow_statistics at the beginning of
847!-- run_control
848    CALL run_control
849!
850!-- Data exchange between coupled models in case that a call has been omitted
851!-- at the end of the previous run of a job chain.
852    IF ( coupling_mode /= 'uncoupled'  .AND.  run_coupled  .AND. .NOT. vnested )  THEN
853!
854!--    In case of model termination initiated by the local model the coupler
855!--    must not be called because this would again cause an MPI hang.
856       DO WHILE ( time_coupling >= dt_coupling  .AND.  terminate_coupled == 0 )
857          CALL surface_coupler
858          time_coupling = time_coupling - dt_coupling
859       ENDDO
860       IF (time_coupling == 0.0_wp  .AND.  time_since_reference_point < dt_coupling )  THEN
861          time_coupling = time_since_reference_point
862       ENDIF
863    ENDIF
864
865    CALL location_message( 'atmosphere (and/or ocean) time-stepping', 'start' )
866
867!
868!-- Start of the time loop
869    DO  WHILE ( simulated_time < end_time  .AND.  .NOT. stop_dt  .AND. .NOT. terminate_run )
870
871       CALL cpu_log( log_point_s(10), 'timesteps', 'start' )
872!
873!--    Vertical nesting: initialize fine grid
874       IF ( vnested ) THEN
875          IF ( .NOT. vnest_init  .AND.  simulated_time >= vnest_start_time )  THEN
876             CALL cpu_log( log_point_s(22), 'vnest_init', 'start' )
877             CALL vnest_init_fine
878             vnest_init = .TRUE.
879             CALL cpu_log( log_point_s(22), 'vnest_init', 'stop' )
880          ENDIF
881       ENDIF
882!
883!--    Determine ug, vg and w_subs in dependence on data from external file
884!--    LSF_DATA
885       IF ( large_scale_forcing .AND. lsf_vert )  THEN
886           CALL ls_forcing_vert ( simulated_time )
887           sums_ls_l = 0.0_wp
888       ENDIF
889
890!
891!--    Set pt_init and q_init to the current profiles taken from
892!--    NUDGING_DATA
893       IF ( nudging )  THEN
894           CALL nudge_ref ( simulated_time )
895!
896!--        Store temperature gradient at the top boundary for possible Neumann
897!--        boundary condition
898           bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
899           bc_q_t_val  = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
900           IF ( air_chemistry )  THEN
901              DO  lsp = 1, nvar
902                 bc_cs_t_val = (  chem_species(lsp)%conc_pr_init(nzt+1)                            &
903                                - chem_species(lsp)%conc_pr_init(nzt) )                            &
904                               / dzu(nzt+1)
905              ENDDO
906           ENDIF
907           IF ( salsa  .AND.  time_since_reference_point >= skip_time_do_salsa )  THEN
908              DO  ib = 1, nbins_aerosol
909                 bc_an_t_val = ( aerosol_number(ib)%init(nzt+1) - aerosol_number(ib)%init(nzt) ) / &
910                               dzu(nzt+1)
911                 DO  ic = 1, ncomponents_mass
912                    icc = ( ic - 1 ) * nbins_aerosol + ib
913                    bc_am_t_val = ( aerosol_mass(icc)%init(nzt+1) - aerosol_mass(icc)%init(nzt) ) /&
914                                  dzu(nzt+1)
915                 ENDDO
916              ENDDO
917              IF ( .NOT. salsa_gases_from_chem )  THEN
918                 DO  ig = 1, ngases_salsa
919                    bc_gt_t_val = ( salsa_gas(ig)%init(nzt+1) - salsa_gas(ig)%init(nzt) ) /        &
920                                  dzu(nzt+1)
921                 ENDDO
922              ENDIF
923           ENDIF
924       ENDIF
925!
926!--    If forcing by larger-scale models is applied, check if new data
927!--    at domain boundaries need to be read.
928       IF ( nesting_offline ) THEN
929          IF ( nest_offl%time(nest_offl%tind_p) <= time_since_reference_point ) &
930               CALL netcdf_data_input_offline_nesting
931       ENDIF
932
933!
934!--    Execute all other module actions routunes
935       CALL module_interface_actions( 'before_timestep' )
936       
937!
938!--    Start of intermediate step loop
939       intermediate_timestep_count = 0
940       DO  WHILE ( intermediate_timestep_count < intermediate_timestep_count_max )
941
942          intermediate_timestep_count = intermediate_timestep_count + 1
943
944!
945!--       Set the steering factors for the prognostic equations which depend
946!--       on the timestep scheme
947          CALL timestep_scheme_steering
948
949!
950!--       Calculate those variables needed in the tendency terms which need
951!--       global communication
952          IF ( .NOT. use_single_reference_value  .AND.  .NOT. use_initial_profile_as_reference )   &
953          THEN
954!
955!--          Horizontally averaged profiles to be used as reference state in
956!--          buoyancy terms (WARNING: only the respective last call of
957!--          calc_mean_profile defines the reference state!)
958             IF ( .NOT. neutral )  THEN
959                CALL calc_mean_profile( pt, 4 )
960                ref_state(:)  = hom(:,1,4,0) ! this is used in the buoyancy term
961             ENDIF
962             IF ( ocean_mode )  THEN
963                CALL calc_mean_profile( rho_ocean, 64 )
964                ref_state(:)  = hom(:,1,64,0)
965             ENDIF
966             IF ( humidity )  THEN
967                CALL calc_mean_profile( vpt, 44 )
968                ref_state(:)  = hom(:,1,44,0)
969             ENDIF
970!
971!--          Assure that ref_state does not become zero at any level
972!--          ( might be the case if a vertical level is completely occupied
973!--            with topography ).
974             ref_state = MERGE( MAXVAL(ref_state), ref_state, ref_state == 0.0_wp )
975          ENDIF
976
977          IF ( .NOT. constant_diffusion )  CALL production_e_init
978          IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  intermediate_timestep_count == 1 )     &
979          THEN
980             CALL ws_statistics
981          ENDIF
982!
983!--       In case of nudging calculate current nudging time scale and horizontal
984!--       means of u, v, pt and q
985          IF ( nudging )  THEN
986             CALL calc_tnudge( simulated_time )
987             CALL calc_mean_profile( u, 1 )
988             CALL calc_mean_profile( v, 2 )
989             CALL calc_mean_profile( pt, 4 )
990             CALL calc_mean_profile( q, 41 )
991          ENDIF
992!
993!--       Execute all other module actions routunes
994          CALL module_interface_actions( 'before_prognostic_equations' )
995!
996!--       Solve the prognostic equations. A fast cache optimized version with
997!--       only one single loop is used in case of Piascek-Williams advection
998!--       scheme. NEC vector machines use a different version, because
999!--       in the other versions a good vectorization is prohibited due to
1000!--       inlining problems.
1001          IF ( loop_optimization == 'cache' )  THEN
1002             CALL prognostic_equations_cache
1003          ELSEIF ( loop_optimization == 'vector' )  THEN
1004             CALL prognostic_equations_vector
1005          ENDIF
1006
1007!
1008!--       Particle transport/physics with the Lagrangian particle model
1009!--       (only once during intermediate steps, because it uses an Euler-step)
1010!--       ### particle model should be moved before prognostic_equations, in order
1011!--       to regard droplet interactions directly
1012          IF ( particle_advection  .AND.  time_since_reference_point >= particle_advection_start   &
1013               .AND.  intermediate_timestep_count == 1 )                                           &
1014          THEN
1015             CALL module_interface_actions( 'after_prognostic_equations' )
1016             first_call_lpm = .FALSE.
1017          ENDIF
1018
1019!
1020!--       Interaction of droplets with temperature and mixing ratio.
1021!--       Droplet condensation and evaporation is calculated within
1022!--       advec_particles.
1023          IF ( cloud_droplets  .AND.  intermediate_timestep_count == intermediate_timestep_count_max ) &
1024          THEN
1025             CALL lpm_interaction_droplets_ptq
1026          ENDIF
1027
1028!
1029!--       Movement of agents in multi agent system
1030          IF ( agents_active  .AND.  time_since_reference_point >= multi_agent_system_start  .AND. &
1031               time_since_reference_point <= multi_agent_system_end  .AND.                         &
1032               intermediate_timestep_count == 1 )                                                  &
1033          THEN
1034             CALL multi_agent_system
1035             first_call_mas = .FALSE.
1036          ENDIF
1037
1038!
1039!--       Exchange of ghost points (lateral boundary conditions)
1040          CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'start' )
1041
1042          CALL exchange_horiz( u_p, nbgp )
1043          CALL exchange_horiz( v_p, nbgp )
1044          CALL exchange_horiz( w_p, nbgp )
1045          CALL exchange_horiz( pt_p, nbgp )
1046          IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e_p, nbgp )
1047          IF ( rans_tke_e  .OR.  wang_kernel  .OR.  collision_turbulence                           &
1048               .OR.  use_sgs_for_particles )  THEN
1049             IF ( rans_tke_e )  THEN
1050                CALL exchange_horiz( diss_p, nbgp )
1051             ELSE
1052                CALL exchange_horiz( diss, nbgp )
1053             ENDIF
1054          ENDIF
1055          IF ( ocean_mode )  THEN
1056             CALL exchange_horiz( sa_p, nbgp )
1057             CALL exchange_horiz( rho_ocean, nbgp )
1058             CALL exchange_horiz( prho, nbgp )
1059          ENDIF
1060          IF ( humidity )  THEN
1061             CALL exchange_horiz( q_p, nbgp )
1062             IF ( bulk_cloud_model .AND. microphysics_morrison )  THEN
1063                CALL exchange_horiz( qc_p, nbgp )
1064                CALL exchange_horiz( nc_p, nbgp )
1065             ENDIF
1066             IF ( bulk_cloud_model .AND. microphysics_seifert )  THEN
1067                CALL exchange_horiz( qr_p, nbgp )
1068                CALL exchange_horiz( nr_p, nbgp )
1069             ENDIF
1070          ENDIF
1071          IF ( cloud_droplets )  THEN
1072             CALL exchange_horiz( ql, nbgp )
1073             CALL exchange_horiz( ql_c, nbgp )
1074             CALL exchange_horiz( ql_v, nbgp )
1075             CALL exchange_horiz( ql_vp, nbgp )
1076          ENDIF
1077          IF ( passive_scalar )  CALL exchange_horiz( s_p, nbgp )
1078          IF ( air_chemistry )  THEN
1079             DO  lsp = 1, nvar
1080                CALL exchange_horiz( chem_species(lsp)%conc_p, nbgp )
1081             ENDDO
1082          ENDIF
1083
1084          IF ( salsa  .AND.  time_since_reference_point >= skip_time_do_salsa )  THEN
1085             DO  ib = 1, nbins_aerosol
1086                CALL exchange_horiz( aerosol_number(ib)%conc_p, nbgp )
1087                DO  ic = 1, ncomponents_mass
1088                   icc = ( ic - 1 ) * nbins_aerosol + ib
1089                   CALL exchange_horiz( aerosol_mass(icc)%conc_p, nbgp )
1090                ENDDO
1091             ENDDO
1092             IF ( .NOT. salsa_gases_from_chem )  THEN
1093                DO  ig = 1, ngases_salsa
1094                   CALL exchange_horiz( salsa_gas(ig)%conc_p, nbgp )
1095                ENDDO
1096             ENDIF
1097          ENDIF
1098          CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
1099
1100!
1101!--       Boundary conditions for the prognostic quantities (except of the
1102!--       velocities at the outflow in case of a non-cyclic lateral wall)
1103          CALL boundary_conds
1104
1105!
1106!--       Boundary conditions for prognostic quantitites of other modules:
1107!--       Here, only decycling is carried out
1108          IF ( air_chemistry )  THEN
1109
1110             DO  lsp = 1, nvar
1111                lsp_usr = 1
1112                DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
1113                   IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN
1114                      CALL chem_boundary_conds( chem_species(lsp)%conc_p,                          &
1115                                                chem_species(lsp)%conc_pr_init )
1116                   ENDIF
1117                   lsp_usr = lsp_usr + 1
1118                ENDDO
1119             ENDDO
1120
1121          ENDIF
1122
1123          IF ( salsa  .AND.  time_since_reference_point >= skip_time_do_salsa )  THEN
1124
1125             DO  ib = 1, nbins_aerosol
1126                CALL salsa_boundary_conds( aerosol_number(ib)%conc_p, aerosol_number(ib)%init )
1127                DO  ic = 1, ncomponents_mass
1128                   icc = ( ic - 1 ) * nbins_aerosol + ib
1129                   CALL salsa_boundary_conds( aerosol_mass(icc)%conc_p, aerosol_mass(icc)%init )
1130                ENDDO
1131             ENDDO
1132             IF ( .NOT. salsa_gases_from_chem )  THEN
1133                DO  ig = 1, ngases_salsa
1134                   CALL salsa_boundary_conds( salsa_gas(ig)%conc_p, salsa_gas(ig)%init )
1135                ENDDO
1136             ENDIF
1137
1138          ENDIF
1139
1140!
1141!--       Incrementing timestep counter
1142          timestep_count = timestep_count + 1
1143
1144          CALL cpu_log( log_point(28), 'swap_timelevel', 'start' )
1145!
1146!--       Set the swap level for all modules
1147          CALL module_interface_swap_timelevel( MOD( timestep_count, 2) )
1148!
1149!--       Set the swap level for steering the pmc data transfer
1150          IF ( nested_run )  CALL pmci_set_swaplevel( MOD( timestep_count, 2) + 1 )  !> @todo: why the +1 ?
1151
1152          CALL cpu_log( log_point(28), 'swap_timelevel', 'stop' )
1153
1154!
1155!--       Vertical nesting: Interpolate fine grid data to the coarse grid
1156          IF ( vnest_init ) THEN
1157             CALL cpu_log( log_point_s(37), 'vnest_anterpolate', 'start' )
1158             CALL vnest_anterpolate
1159             CALL cpu_log( log_point_s(37), 'vnest_anterpolate', 'stop' )
1160          ENDIF
1161
1162          IF ( nested_run )  THEN
1163
1164             CALL cpu_log( log_point(60), 'nesting', 'start' )
1165!
1166!--          Domain nesting. The data transfer subroutines pmci_parent_datatrans
1167!--          and pmci_child_datatrans are called inside the wrapper
1168!--          subroutine pmci_datatrans according to the control parameters
1169!--          nesting_mode and nesting_datatransfer_mode.
1170!--          TO_DO: why is nesting_mode given as a parameter here?
1171             CALL pmci_datatrans( nesting_mode )
1172
1173             IF ( TRIM( nesting_mode ) == 'two-way' .OR.  nesting_mode == 'vertical' )  THEN
1174
1175                CALL cpu_log( log_point_s(92), 'exchange-horiz-nest', 'start' )
1176!
1177!--             Exchange_horiz is needed for all parent-domains after the
1178!--             anterpolation
1179                CALL exchange_horiz( u, nbgp )
1180                CALL exchange_horiz( v, nbgp )
1181                CALL exchange_horiz( w, nbgp )
1182                IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
1183
1184                IF ( humidity )  THEN
1185
1186                   CALL exchange_horiz( q, nbgp )
1187
1188                   IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
1189                       CALL exchange_horiz( qc, nbgp )
1190                       CALL exchange_horiz( nc, nbgp )
1191                   ENDIF
1192                   IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
1193                       CALL exchange_horiz( qr, nbgp )
1194                       CALL exchange_horiz( nr, nbgp )
1195                   ENDIF
1196
1197                ENDIF
1198
1199                IF ( passive_scalar )  CALL exchange_horiz( s, nbgp ) 
1200
1201                IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e, nbgp )
1202
1203                IF ( .NOT. constant_diffusion  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
1204                   CALL exchange_horiz( diss, nbgp )
1205                ENDIF
1206
1207                IF ( air_chemistry )  THEN
1208                   DO  n = 1, nvar
1209                      CALL exchange_horiz( chem_species(n)%conc, nbgp ) 
1210                   ENDDO
1211                ENDIF
1212
1213                IF ( salsa  .AND. time_since_reference_point >= skip_time_do_salsa )  THEN
1214                   DO  ib = 1, nbins_aerosol
1215                      CALL exchange_horiz( aerosol_number(ib)%conc, nbgp )
1216                      DO  ic = 1, ncomponents_mass
1217                         icc = ( ic - 1 ) * nbins_aerosol + ib
1218                         CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp )
1219                      ENDDO
1220                   ENDDO
1221                   IF ( .NOT. salsa_gases_from_chem )  THEN
1222                      DO  ig = 1, ngases_salsa
1223                         CALL exchange_horiz( salsa_gas(ig)%conc, nbgp )
1224                      ENDDO
1225                   ENDIF
1226                ENDIF
1227                CALL cpu_log( log_point_s(92), 'exchange-horiz-nest', 'stop' )
1228
1229             ENDIF
1230
1231!
1232!--          Set boundary conditions again after interpolation and anterpolation.
1233             CALL pmci_boundary_conds
1234
1235!
1236!--          Set chemistry boundary conditions (decycling)
1237             IF ( air_chemistry )  THEN
1238                DO  lsp = 1, nvar
1239                   lsp_usr = 1
1240                   DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
1241                      IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN
1242                         CALL chem_boundary_conds( chem_species(lsp)%conc,                         &
1243                                                   chem_species(lsp)%conc_pr_init )
1244                      ENDIF
1245                      lsp_usr = lsp_usr + 1
1246                   ENDDO
1247                ENDDO
1248             ENDIF
1249
1250!
1251!--          Set SALSA boundary conditions (decycling)
1252             IF ( salsa  .AND. time_since_reference_point >= skip_time_do_salsa )  THEN
1253                DO  ib = 1, nbins_aerosol
1254                   CALL salsa_boundary_conds( aerosol_number(ib)%conc, aerosol_number(ib)%init )
1255                   DO  ic = 1, ncomponents_mass
1256                      icc = ( ic - 1 ) * nbins_aerosol + ib
1257                      CALL salsa_boundary_conds( aerosol_mass(icc)%conc, aerosol_mass(icc)%init )
1258                   ENDDO
1259                ENDDO
1260                IF ( .NOT. salsa_gases_from_chem )  THEN
1261                   DO  ig = 1, ngases_salsa
1262                      CALL salsa_boundary_conds( salsa_gas(ig)%conc, salsa_gas(ig)%init )
1263                   ENDDO
1264                ENDIF
1265             ENDIF
1266
1267             CALL cpu_log( log_point(60), 'nesting', 'stop' )
1268
1269          ENDIF
1270
1271!
1272!--       Temperature offset must be imposed at cyclic boundaries in x-direction
1273!--       when a sloping surface is used
1274          IF ( sloping_surface )  THEN
1275             IF ( nxl ==  0 )  pt(:,:,nxlg:nxl-1) = pt(:,:,nxlg:nxl-1) - pt_slope_offset
1276             IF ( nxr == nx )  pt(:,:,nxr+1:nxrg) = pt(:,:,nxr+1:nxrg) + pt_slope_offset
1277          ENDIF
1278
1279!
1280!--       Impose a turbulent inflow using the recycling method
1281          IF ( turbulent_inflow )  CALL inflow_turbulence
1282
1283!
1284!--       Set values at outflow boundary using the special outflow condition
1285          IF ( turbulent_outflow )  CALL outflow_turbulence
1286
1287!
1288!--       Impose a random perturbation on the horizontal velocity field
1289          IF ( create_disturbances  .AND.  ( call_psolver_at_all_substeps  .AND.                   &
1290               intermediate_timestep_count == intermediate_timestep_count_max )                    &
1291               .OR. ( .NOT. call_psolver_at_all_substeps  .AND.  intermediate_timestep_count == 1 ) ) &
1292          THEN
1293             time_disturb = time_disturb + dt_3d
1294             IF ( time_disturb >= dt_disturb )  THEN
1295                IF ( disturbance_energy_limit /= 0.0_wp  .AND.                                     &
1296                     hom(nzb+5,1,pr_palm,0) < disturbance_energy_limit )  THEN
1297                   CALL disturb_field( 'u', tend, u )
1298                   CALL disturb_field( 'v', tend, v )
1299                ELSEIF ( ( .NOT. bc_lr_cyc  .OR.  .NOT. bc_ns_cyc )                                &
1300                         .AND. .NOT. child_domain  .AND.  .NOT.  nesting_offline )                 &
1301                THEN
1302!
1303!--                Runs with a non-cyclic lateral wall need perturbations
1304!--                near the inflow throughout the whole simulation
1305                   dist_range = 1
1306                   CALL disturb_field( 'u', tend, u )
1307                   CALL disturb_field( 'v', tend, v )
1308                   dist_range = 0
1309                ENDIF
1310                time_disturb = time_disturb - dt_disturb
1311             ENDIF
1312          ENDIF
1313
1314!
1315!--       Map forcing data derived from larger scale model onto domain
1316!--       boundaries.
1317          IF ( nesting_offline  .AND.  intermediate_timestep_count ==                              &
1318                                       intermediate_timestep_count_max  )                          &
1319             CALL nesting_offl_bc
1320!
1321!--       Impose a turbulent inflow using synthetic generated turbulence.
1322          IF ( use_syn_turb_gen  .AND.                                                             &
1323               intermediate_timestep_count == intermediate_timestep_count_max )  THEN
1324             CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
1325             CALL stg_main
1326             CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
1327          ENDIF
1328!
1329!--       Ensure mass conservation. This need to be done after imposing
1330!--       synthetic turbulence and top boundary condition for pressure is set to
1331!--       Neumann conditions.
1332!--       Is this also required in case of Dirichlet?
1333          IF ( nesting_offline )  CALL nesting_offl_mass_conservation
1334!
1335!--       Reduce the velocity divergence via the equation for perturbation
1336!--       pressure.
1337          IF ( intermediate_timestep_count == 1  .OR. &
1338                call_psolver_at_all_substeps )  THEN
1339
1340             IF (  vnest_init ) THEN
1341!
1342!--             Compute pressure in the CG, interpolate top boundary conditions
1343!--             to the FG and then compute pressure in the FG
1344                IF ( coupling_mode == 'vnested_crse' )  CALL pres
1345
1346                CALL cpu_log( log_point_s(30), 'vnest_bc', 'start' )
1347                CALL vnest_boundary_conds
1348                CALL cpu_log( log_point_s(30), 'vnest_bc', 'stop' )
1349 
1350                IF ( coupling_mode == 'vnested_fine' )  CALL pres
1351
1352!--             Anterpolate TKE, satisfy Germano Identity
1353                CALL cpu_log( log_point_s(28), 'vnest_anter_e', 'start' )
1354                CALL vnest_anterpolate_e
1355                CALL cpu_log( log_point_s(28), 'vnest_anter_e', 'stop' )
1356
1357             ELSE
1358!               
1359!--             Mass (volume) flux correction to ensure global mass conservation for child domains.
1360                IF ( child_domain )  THEN
1361                   IF ( nesting_mode == 'vertical' )  THEN
1362                      CALL pmci_ensure_nest_mass_conservation_vertical
1363                   ELSE
1364                      CALL pmci_ensure_nest_mass_conservation
1365                   ENDIF
1366                ENDIF
1367               
1368                CALL pres
1369
1370             ENDIF
1371
1372          ENDIF
1373
1374!
1375!--       If required, compute liquid water content
1376          IF ( bulk_cloud_model )  THEN
1377             CALL calc_liquid_water_content
1378          ENDIF
1379!
1380!--       If required, compute virtual potential temperature
1381          IF ( humidity )  THEN
1382             CALL compute_vpt
1383          ENDIF
1384
1385!
1386!--       Compute the diffusion quantities
1387          IF ( .NOT. constant_diffusion )  THEN
1388
1389!
1390!--          Determine surface fluxes shf and qsws and surface values
1391!--          pt_surface and q_surface in dependence on data from external
1392!--          file LSF_DATA respectively
1393             IF ( ( large_scale_forcing .AND. lsf_surf ) .AND.                                     &
1394                 intermediate_timestep_count == intermediate_timestep_count_max )                  &
1395             THEN
1396                CALL ls_forcing_surf( simulated_time )
1397             ENDIF
1398
1399!
1400!--          First the vertical (and horizontal) fluxes in the surface
1401!--          (constant flux) layer are computed
1402             IF ( constant_flux_layer )  THEN
1403                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'start' )
1404                CALL surface_layer_fluxes
1405                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'stop' )
1406             ENDIF
1407!
1408!--          If required, solve the energy balance for the surface and run soil
1409!--          model. Call for horizontal as well as vertical surfaces
1410             IF ( land_surface .AND. time_since_reference_point >= skip_time_do_lsm)  THEN
1411
1412                CALL cpu_log( log_point(54), 'land_surface', 'start' )
1413!
1414!--             Call for horizontal upward-facing surfaces
1415                CALL lsm_energy_balance( .TRUE., -1 )
1416                CALL lsm_soil_model( .TRUE., -1, .TRUE. )
1417!
1418!--             Call for northward-facing surfaces
1419                CALL lsm_energy_balance( .FALSE., 0 )
1420                CALL lsm_soil_model( .FALSE., 0, .TRUE. )
1421!
1422!--             Call for southward-facing surfaces
1423                CALL lsm_energy_balance( .FALSE., 1 )
1424                CALL lsm_soil_model( .FALSE., 1, .TRUE. )
1425!
1426!--             Call for eastward-facing surfaces
1427                CALL lsm_energy_balance( .FALSE., 2 )
1428                CALL lsm_soil_model( .FALSE., 2, .TRUE. )
1429!
1430!--             Call for westward-facing surfaces
1431                CALL lsm_energy_balance( .FALSE., 3 )
1432                CALL lsm_soil_model( .FALSE., 3, .TRUE. )
1433               
1434!
1435!--             At the end, set boundary conditons for potential temperature
1436!--             and humidity after running the land-surface model. This
1437!--             might be important for the nesting, where arrays are transfered.
1438                CALL lsm_boundary_condition
1439
1440               
1441                CALL cpu_log( log_point(54), 'land_surface', 'stop' )
1442             ENDIF
1443!
1444!--          If required, solve the energy balance for urban surfaces and run
1445!--          the material heat model
1446             IF (urban_surface) THEN
1447                CALL cpu_log( log_point(74), 'urban_surface', 'start' )
1448               
1449                CALL usm_surface_energy_balance( .FALSE. )
1450                IF ( usm_material_model )  THEN
1451                   CALL usm_green_heat_model
1452                   CALL usm_material_heat_model ( .FALSE. )
1453                ENDIF
1454
1455!
1456!--             At the end, set boundary conditons for potential temperature
1457!--             and humidity after running the urban-surface model. This
1458!--             might be important for the nesting, where arrays are transfered.
1459                CALL usm_boundary_condition
1460
1461                CALL cpu_log( log_point(74), 'urban_surface', 'stop' )
1462             ENDIF
1463!
1464!--          Compute the diffusion coefficients
1465             CALL cpu_log( log_point(17), 'diffusivities', 'start' )
1466             IF ( .NOT. humidity ) THEN
1467                IF ( ocean_mode )  THEN
1468                   CALL tcm_diffusivities( prho, prho_reference )
1469                ELSE
1470                   CALL tcm_diffusivities( pt, pt_reference )
1471                ENDIF
1472             ELSE
1473                CALL tcm_diffusivities( vpt, pt_reference )
1474             ENDIF
1475             CALL cpu_log( log_point(17), 'diffusivities', 'stop' )
1476!
1477!--          Vertical nesting: set fine grid eddy viscosity top boundary condition
1478             IF ( vnest_init )  CALL vnest_boundary_conds_khkm
1479
1480          ENDIF
1481
1482!
1483!--       If required, calculate radiative fluxes and heating rates
1484          IF ( radiation  .AND.  intermediate_timestep_count == intermediate_timestep_count_max    &
1485               .AND. time_since_reference_point > skip_time_do_radiation )  THEN
1486
1487               time_radiation = time_radiation + dt_3d
1488
1489             IF ( time_radiation >= dt_radiation  .OR.  force_radiation_call )  THEN
1490
1491                CALL cpu_log( log_point(50), 'radiation', 'start' )
1492
1493                IF ( .NOT. force_radiation_call )  THEN
1494                   time_radiation = time_radiation - dt_radiation
1495                ENDIF
1496
1497!
1498!--             Adjust the current time to the time step of the radiation model.
1499!--             Needed since radiation is pre-calculated and stored only on apparent
1500!--             solar positions
1501                time_since_reference_point_save = time_since_reference_point
1502                time_since_reference_point = REAL( FLOOR( time_since_reference_point /             &
1503                                                          dt_radiation), wp ) * dt_radiation
1504
1505                CALL radiation_control
1506
1507                IF ( ( urban_surface  .OR.  land_surface )  .AND.  radiation_interactions )  THEN
1508                   CALL cpu_log( log_point_s(46), 'radiation_interaction', 'start' )
1509                   CALL radiation_interaction
1510                   CALL cpu_log( log_point_s(46), 'radiation_interaction', 'stop' )
1511                ENDIF
1512   
1513!
1514!--             Return the current time to its original value
1515                time_since_reference_point = time_since_reference_point_save
1516
1517                CALL cpu_log( log_point(50), 'radiation', 'stop' )
1518
1519             ENDIF
1520          ENDIF
1521
1522       ENDDO   ! Intermediate step loop
1523
1524!
1525!--    Will be used at some point by flow_statistics.
1526       !$ACC UPDATE &
1527       !$ACC HOST(sums_l_l(nzb:nzt+1,0:statistic_regions,0)) &
1528       !$ACC HOST(sums_us2_ws_l(nzb:nzt+1,0)) &
1529       !$ACC HOST(sums_wsus_ws_l(nzb:nzt+1,0)) &
1530       !$ACC HOST(sums_vs2_ws_l(nzb:nzt+1,0)) &
1531       !$ACC HOST(sums_wsvs_ws_l(nzb:nzt+1,0)) &
1532       !$ACC HOST(sums_ws2_ws_l(nzb:nzt+1,0)) &
1533       !$ACC HOST(sums_wspts_ws_l(nzb:nzt+1,0)) &
1534       !$ACC HOST(sums_wssas_ws_l(nzb:nzt+1,0)) &
1535       !$ACC HOST(sums_wsqs_ws_l(nzb:nzt+1,0)) &
1536       !$ACC HOST(sums_wsqcs_ws_l(nzb:nzt+1,0)) &
1537       !$ACC HOST(sums_wsqrs_ws_l(nzb:nzt+1,0)) &
1538       !$ACC HOST(sums_wsncs_ws_l(nzb:nzt+1,0)) &
1539       !$ACC HOST(sums_wsnrs_ws_l(nzb:nzt+1,0)) &
1540       !$ACC HOST(sums_wsss_ws_l(nzb:nzt+1,0)) &
1541       !$ACC HOST(sums_salsa_ws_l(nzb:nzt+1,0))
1542
1543!
1544!--    If required, consider chemical emissions
1545       IF ( air_chemistry  .AND.  emissions_anthropogenic )  THEN
1546!
1547!--       Update the time --> kanani: revise location of this CALL
1548          CALL calc_date_and_time
1549!
1550!--       Call emission routine only once an hour
1551          IF (hour_of_year  .GT.  hour_call_emis )  THEN
1552             CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
1553             hour_call_emis = hour_of_year
1554          ENDIF
1555       ENDIF
1556!
1557!--    If required, consider aerosol emissions for the salsa model
1558       IF ( salsa )  THEN
1559!
1560!--       Call emission routine to update emissions if needed
1561          CALL salsa_emission_update
1562
1563       ENDIF
1564!
1565!--    If required, calculate indoor temperature, waste heat, heat flux
1566!--    through wall, etc.
1567!--    dt_indoor steers the frequency of the indoor model calculations.
1568!--    Note, at first timestep indoor model is called, in order to provide
1569!--    a waste heat flux.
1570       IF ( indoor_model )  THEN
1571
1572          time_indoor = time_indoor + dt_3d
1573
1574          IF ( time_indoor >= dt_indoor  .OR.  current_timestep_number == 0 )  THEN
1575
1576             time_indoor = time_indoor - dt_indoor
1577
1578             CALL cpu_log( log_point(76), 'indoor_model', 'start' )
1579             CALL im_main_heatcool
1580             CALL cpu_log( log_point(76), 'indoor_model', 'stop' )
1581
1582          ENDIF
1583       ENDIF
1584!
1585!--    Increase simulation time and output times
1586       nr_timesteps_this_run      = nr_timesteps_this_run + 1
1587       current_timestep_number    = current_timestep_number + 1
1588       simulated_time             = simulated_time   + dt_3d
1589       time_since_reference_point = simulated_time - coupling_start_time
1590       simulated_time_chr         = time_to_string( time_since_reference_point )
1591
1592
1593
1594
1595       IF ( time_since_reference_point >= skip_time_data_output_av )  THEN
1596          time_do_av         = time_do_av       + dt_3d
1597       ENDIF
1598       IF ( time_since_reference_point >= skip_time_do2d_xy )  THEN
1599          time_do2d_xy       = time_do2d_xy     + dt_3d
1600       ENDIF
1601       IF ( time_since_reference_point >= skip_time_do2d_xz )  THEN
1602          time_do2d_xz       = time_do2d_xz     + dt_3d
1603       ENDIF
1604       IF ( time_since_reference_point >= skip_time_do2d_yz )  THEN
1605          time_do2d_yz       = time_do2d_yz     + dt_3d
1606       ENDIF
1607       IF ( time_since_reference_point >= skip_time_do3d    )  THEN
1608          time_do3d          = time_do3d        + dt_3d
1609       ENDIF
1610       DO  mid = 1, masks
1611          IF ( time_since_reference_point >= skip_time_domask(mid) )  THEN
1612             time_domask(mid)= time_domask(mid) + dt_3d
1613          ENDIF
1614       ENDDO
1615       IF ( time_since_reference_point >= skip_time_dosp )  THEN
1616          time_dosp       = time_dosp        + dt_3d
1617       ENDIF
1618       time_dots          = time_dots        + dt_3d
1619       IF ( .NOT. first_call_lpm )  THEN
1620          time_dopts      = time_dopts       + dt_3d
1621       ENDIF
1622       IF ( time_since_reference_point >= skip_time_dopr )  THEN
1623          time_dopr       = time_dopr        + dt_3d
1624       ENDIF
1625       time_dopr_listing  = time_dopr_listing + dt_3d
1626       time_run_control   = time_run_control + dt_3d
1627!
1628!--    Increment time-counter for surface output
1629       IF ( surface_output )  THEN
1630          IF ( time_since_reference_point >= skip_time_dosurf )  THEN
1631             time_dosurf    = time_dosurf + dt_3d
1632          ENDIF
1633          IF ( time_since_reference_point >= skip_time_dosurf_av )  THEN
1634             time_dosurf_av = time_dosurf_av + dt_3d
1635          ENDIF
1636       ENDIF
1637!
1638!--    Increment time-counter for virtual measurements
1639       IF ( virtual_measurement  .AND.  vm_time_start <= time_since_reference_point )  THEN
1640          time_virtual_measurement = time_virtual_measurement + dt_3d
1641       ENDIF
1642!
1643!--    In case of synthetic turbulence generation and parametrized turbulence
1644!--    information, update the time counter and if required, adjust the
1645!--    STG to new atmospheric conditions.
1646       IF ( use_syn_turb_gen  )  THEN
1647          IF ( parametrize_inflow_turbulence )  THEN
1648             time_stg_adjust = time_stg_adjust + dt_3d
1649             IF ( time_stg_adjust >= dt_stg_adjust )  THEN
1650                CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
1651                CALL stg_adjust
1652                CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
1653             ENDIF
1654          ENDIF
1655          time_stg_call = time_stg_call + dt_3d
1656       ENDIF
1657
1658!
1659!--    Data exchange between coupled models
1660       IF ( coupling_mode /= 'uncoupled'  .AND.  run_coupled  .AND.  .NOT. vnested )  THEN
1661          time_coupling = time_coupling + dt_3d
1662
1663!
1664!--       In case of model termination initiated by the local model
1665!--       (terminate_coupled > 0), the coupler must be skipped because it would
1666!--       cause an MPI intercomminucation hang.
1667!--       If necessary, the coupler will be called at the beginning of the
1668!--       next restart run.
1669          DO WHILE ( time_coupling >= dt_coupling  .AND.  terminate_coupled == 0 )
1670             CALL surface_coupler
1671             time_coupling = time_coupling - dt_coupling
1672          ENDDO
1673       ENDIF
1674
1675!
1676!--    Biometeorology calculation of stationary thermal indices
1677!--    Todo (kanani): biometeorology needs own time_... treatment.
1678!--                   It might be that time_do2d_xy differs from time_do3d,
1679!--                   and then we might get trouble with the biomet output,
1680!--                   because we can have 2d and/or 3d biomet output!!
1681       IF ( biometeorology                                                                         &
1682            .AND. ( ( time_do3d >= dt_do3d  .AND.  time_since_reference_point >= skip_time_do3d )  &
1683                  .OR.                                                                             &
1684            ( time_do2d_xy >= dt_do2d_xy  .AND.  time_since_reference_point >= skip_time_do2d_xy ) &
1685                    ) )  THEN
1686!
1687!--       If required, do thermal comfort calculations
1688          IF ( thermal_comfort )  THEN
1689             CALL bio_calculate_thermal_index_maps ( .FALSE. )
1690          ENDIF
1691!
1692!--       If required, do UV exposure calculations
1693          IF ( uv_exposure )  THEN
1694             CALL uvem_calc_exposure
1695          ENDIF
1696       ENDIF
1697
1698!
1699!--    Execute alle other module actions routunes
1700       CALL module_interface_actions( 'after_integration' )
1701
1702!
1703!--    If Galilei transformation is used, determine the distance that the
1704!--    model has moved so far
1705       IF ( galilei_transformation )  THEN
1706          advected_distance_x = advected_distance_x + u_gtrans * dt_3d
1707          advected_distance_y = advected_distance_y + v_gtrans * dt_3d
1708       ENDIF
1709
1710!
1711!--    Check, if restart is necessary (because cpu-time is expiring or
1712!--    because it is forced by user) and set stop flag
1713!--    This call is skipped if the remote model has already initiated a restart.
1714       IF ( .NOT. terminate_run )  CALL check_for_restart
1715
1716!
1717!--    Carry out statistical analysis and output at the requested output times.
1718!--    The MOD function is used for calculating the output time counters (like
1719!--    time_dopr) in order to regard a possible decrease of the output time
1720!--    interval in case of restart runs
1721
1722!
1723!--    Set a flag indicating that so far no statistics have been created
1724!--    for this time step
1725       flow_statistics_called = .FALSE.
1726
1727!
1728!--    If required, call flow_statistics for averaging in time
1729       IF ( averaging_interval_pr /= 0.0_wp  .AND.                                                 &
1730            ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.                                &
1731            time_since_reference_point >= skip_time_dopr )  THEN
1732          time_dopr_av = time_dopr_av + dt_3d
1733          IF ( time_dopr_av >= dt_averaging_input_pr )  THEN
1734             do_sum = .TRUE.
1735             time_dopr_av = MOD( time_dopr_av, MAX( dt_averaging_input_pr, dt_3d ) )
1736          ENDIF
1737       ENDIF
1738       IF ( do_sum )  CALL flow_statistics
1739
1740!
1741!--    Sum-up 3d-arrays for later output of time-averaged 2d/3d/masked data
1742       IF ( averaging_interval /= 0.0_wp  .AND.                                                    &
1743            ( dt_data_output_av - time_do_av ) <= averaging_interval  .AND.                        &
1744            time_since_reference_point >= skip_time_data_output_av )                               &
1745       THEN
1746          time_do_sla = time_do_sla + dt_3d
1747          IF ( time_do_sla >= dt_averaging_input )  THEN
1748             IF ( current_timestep_number > timestep_number_at_prev_calc )                         &
1749                CALL doq_calculate
1750
1751             CALL sum_up_3d_data
1752             average_count_3d = average_count_3d + 1
1753             time_do_sla = MOD( time_do_sla, MAX( dt_averaging_input, dt_3d ) )
1754          ENDIF
1755       ENDIF
1756!
1757!--    Average surface data
1758       IF ( surface_output )  THEN
1759          IF ( averaging_interval_surf /= 0.0_wp                                                   &
1760                .AND.  ( dt_dosurf_av - time_dosurf_av ) <= averaging_interval_surf                &
1761                .AND.  time_since_reference_point >= skip_time_dosurf_av )  THEN
1762             IF ( time_dosurf_av >= dt_averaging_input )  THEN       
1763                CALL surface_data_output_averaging
1764                average_count_surf = average_count_surf + 1
1765             ENDIF
1766          ENDIF
1767       ENDIF
1768
1769!
1770!--    Calculate spectra for time averaging
1771       IF ( averaging_interval_sp /= 0.0_wp  .AND. ( dt_dosp - time_dosp ) <= averaging_interval_sp&
1772            .AND.  time_since_reference_point >= skip_time_dosp )  THEN
1773          time_dosp_av = time_dosp_av + dt_3d
1774          IF ( time_dosp_av >= dt_averaging_input_pr )  THEN
1775             CALL calc_spectra
1776             time_dosp_av = MOD( time_dosp_av, MAX( dt_averaging_input_pr, dt_3d ) )
1777          ENDIF
1778       ENDIF
1779
1780!
1781!--    Call flight module and output data
1782       IF ( virtual_flight )  THEN
1783          CALL flight_measurement
1784          CALL data_output_flight
1785       ENDIF
1786!
1787!--    Take virtual measurements
1788       IF ( virtual_measurement  .AND.  time_virtual_measurement >= dt_virtual_measurement         &
1789                                 .AND.  vm_time_start <= time_since_reference_point )  THEN
1790          CALL vm_sampling
1791          CALL vm_data_output
1792          time_virtual_measurement = MOD(      time_virtual_measurement,                           &
1793                                          MAX( dt_virtual_measurement, dt_3d ) )
1794       ENDIF
1795!
1796!--    Profile output (ASCII) on file
1797       IF ( time_dopr_listing >= dt_dopr_listing )  THEN
1798          CALL print_1d
1799          time_dopr_listing = MOD( time_dopr_listing, MAX( dt_dopr_listing, dt_3d ) )
1800       ENDIF
1801
1802!
1803!--    Graphic output for PROFIL
1804       IF ( time_dopr >= dt_dopr  .AND.  time_since_reference_point >= skip_time_dopr )  THEN
1805          IF ( dopr_n /= 0 )  CALL data_output_profiles
1806          time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
1807          time_dopr_av = 0.0_wp    ! due to averaging (see above)
1808       ENDIF
1809
1810!
1811!--    Graphic output for time series
1812       IF ( time_dots >= dt_dots )  THEN
1813          CALL data_output_tseries
1814          time_dots = MOD( time_dots, MAX( dt_dots, dt_3d ) )
1815       ENDIF
1816
1817!
1818!--    Output of spectra (formatted for use with PROFIL), in case of no
1819!--    time averaging, spectra has to be calculated before
1820       IF ( time_dosp >= dt_dosp  .AND.  time_since_reference_point >= skip_time_dosp )  THEN
1821          IF ( average_count_sp == 0 )  CALL calc_spectra
1822          CALL data_output_spectra
1823          time_dosp = MOD( time_dosp, MAX( dt_dosp, dt_3d ) )
1824       ENDIF
1825
1826!
1827!--    2d-data output (cross-sections)
1828       IF ( time_do2d_xy >= dt_do2d_xy  .AND.  time_since_reference_point >= skip_time_do2d_xy )  THEN
1829          IF ( current_timestep_number > timestep_number_at_prev_calc )                            &
1830             CALL doq_calculate
1831
1832          CALL data_output_2d( 'xy', 0 )
1833          time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
1834       ENDIF
1835       IF ( time_do2d_xz >= dt_do2d_xz  .AND.  time_since_reference_point >= skip_time_do2d_xz )  THEN
1836          IF ( current_timestep_number > timestep_number_at_prev_calc )                            &
1837
1838             CALL doq_calculate
1839          CALL data_output_2d( 'xz', 0 )
1840          time_do2d_xz = MOD( time_do2d_xz, MAX( dt_do2d_xz, dt_3d ) )
1841       ENDIF
1842       IF ( time_do2d_yz >= dt_do2d_yz  .AND.  time_since_reference_point >= skip_time_do2d_yz )  THEN
1843          IF ( current_timestep_number > timestep_number_at_prev_calc )                            &
1844             CALL doq_calculate
1845
1846          CALL data_output_2d( 'yz', 0 )
1847          time_do2d_yz = MOD( time_do2d_yz, MAX( dt_do2d_yz, dt_3d ) )
1848       ENDIF
1849
1850!
1851!--    3d-data output (volume data)
1852       IF ( time_do3d >= dt_do3d  .AND.  time_since_reference_point >= skip_time_do3d )  THEN
1853          IF ( current_timestep_number > timestep_number_at_prev_calc )                            &
1854             CALL doq_calculate
1855
1856          CALL data_output_3d( 0 )
1857          time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
1858       ENDIF
1859
1860!
1861!--    Masked data output
1862       DO  mid = 1, masks
1863          IF ( time_domask(mid) >= dt_domask(mid)                                                  &
1864               .AND.  time_since_reference_point >= skip_time_domask(mid) )  THEN
1865             IF ( current_timestep_number > timestep_number_at_prev_calc )                         &
1866                CALL doq_calculate
1867
1868             CALL data_output_mask( 0 )
1869             time_domask(mid) = MOD( time_domask(mid), MAX( dt_domask(mid), dt_3d ) )
1870          ENDIF
1871       ENDDO
1872
1873!
1874!--    Output of time-averaged 2d/3d/masked data
1875       IF ( time_do_av >= dt_data_output_av                                                        &
1876            .AND.  time_since_reference_point >= skip_time_data_output_av )  THEN
1877          CALL average_3d_data
1878!
1879!--       Udate thermal comfort indices based on updated averaged input
1880          IF ( biometeorology  .AND.  thermal_comfort )  THEN
1881             CALL bio_calculate_thermal_index_maps ( .TRUE. )
1882          ENDIF
1883          CALL data_output_2d( 'xy', 1 )
1884          CALL data_output_2d( 'xz', 1 )
1885          CALL data_output_2d( 'yz', 1 )
1886          CALL data_output_3d( 1 )
1887          DO  mid = 1, masks
1888             CALL data_output_mask( 1 )
1889          ENDDO
1890          time_do_av = MOD( time_do_av, MAX( dt_data_output_av, dt_3d ) )
1891       ENDIF
1892!
1893!--    Output of surface data, instantaneous and averaged data
1894       IF ( surface_output )  THEN
1895          IF ( time_dosurf >= dt_dosurf  .AND.  time_since_reference_point >= skip_time_dosurf )  THEN
1896             CALL surface_data_output( 0 )
1897             time_dosurf = MOD( time_dosurf, MAX( dt_dosurf, dt_3d ) )
1898          ENDIF
1899          IF ( time_dosurf_av >= dt_dosurf_av  .AND.  time_since_reference_point >= skip_time_dosurf_av )  THEN
1900             CALL surface_data_output( 1 )
1901             time_dosurf_av = MOD( time_dosurf_av, MAX( dt_dosurf_av, dt_3d ) )
1902          ENDIF
1903       ENDIF
1904
1905!
1906!--    Output of particle time series
1907       IF ( particle_advection )  THEN
1908          IF ( time_dopts >= dt_dopts  .OR.                                                        &
1909               ( time_since_reference_point >= particle_advection_start  .AND.                     &
1910                 first_call_lpm ) )  THEN
1911             CALL lpm_data_output_ptseries
1912             time_dopts = MOD( time_dopts, MAX( dt_dopts, dt_3d ) )
1913          ENDIF
1914       ENDIF
1915
1916!
1917!--    If required, set the heat flux for the next time step to a random value
1918       IF ( constant_heatflux  .AND.  random_heatflux )  THEN
1919          IF ( surf_def_h(0)%ns >= 1 )  THEN
1920             CALL cpu_log( log_point(23), 'disturb_heatflux', 'start' )
1921             CALL disturb_heatflux( surf_def_h(0) )
1922             CALL cpu_log( log_point(23), 'disturb_heatflux', 'stop' )
1923          ENDIF
1924          IF ( surf_lsm_h%ns    >= 1 )  THEN
1925             CALL cpu_log( log_point(23), 'disturb_heatflux', 'start' )
1926             CALL disturb_heatflux( surf_lsm_h    )
1927             CALL cpu_log( log_point(23), 'disturb_heatflux', 'stop' )
1928          ENDIF
1929          IF ( surf_usm_h%ns    >= 1 )  THEN
1930             CALL cpu_log( log_point(23), 'disturb_heatflux', 'start' )
1931             CALL disturb_heatflux( surf_usm_h    )
1932             CALL cpu_log( log_point(23), 'disturb_heatflux', 'stop' )
1933          ENDIF
1934       ENDIF
1935
1936!
1937!--    Execute alle other module actions routunes
1938       CALL module_interface_actions( 'after_timestep' )
1939
1940!
1941!--    Determine size of next time step. Save timestep dt_3d because it is
1942!--    newly calculated in routine timestep, but required further below for
1943!--    steering the run control output interval
1944       dt_3d_old = dt_3d
1945       CALL timestep
1946
1947!
1948!--    Synchronize the timestep in case of nested run.
1949       IF ( nested_run )  THEN
1950!
1951!--       Synchronize by unifying the time step.
1952!--       Global minimum of all time-steps is used for all.
1953          CALL pmci_synchronize
1954       ENDIF
1955
1956!
1957!--    Computation and output of run control parameters.
1958!--    This is also done whenever perturbations have been imposed
1959       IF ( time_run_control >= dt_run_control  .OR.                                               &
1960            timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created )                           &
1961       THEN
1962          CALL run_control
1963          IF ( time_run_control >= dt_run_control )  THEN
1964             time_run_control = MOD( time_run_control, MAX( dt_run_control, dt_3d_old ) )
1965          ENDIF
1966       ENDIF
1967
1968!
1969!--    Output elapsed simulated time in form of a progress bar on stdout
1970       IF ( myid == 0 )  CALL output_progress_bar
1971
1972       CALL cpu_log( log_point_s(10), 'timesteps', 'stop' )
1973
1974
1975    ENDDO   ! time loop
1976
1977#if defined( _OPENACC )
1978    CALL exit_surface_arrays
1979#endif
1980!$ACC END DATA
1981!$ACC END DATA
1982!$ACC END DATA
1983!$ACC END DATA
1984!$ACC END DATA
1985!$ACC END DATA
1986!$ACC END DATA
1987
1988!
1989!-- Vertical nesting: Deallocate variables initialized for vertical nesting   
1990    IF ( vnest_init )  CALL vnest_deallocate
1991
1992    IF ( myid == 0 )  CALL finish_progress_bar
1993
1994    CALL location_message( 'atmosphere (and/or ocean) time-stepping', 'finished' )
1995
1996 END SUBROUTINE time_integration
Note: See TracBrowser for help on using the repository browser.