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

Last change on this file since 3994 was 3994, checked in by suehring, 6 years ago

new module for diagnostic output quantities added + output of turbulence intensity

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