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

Last change on this file since 4420 was 4420, checked in by maronga, 4 years ago

added steering for NetCDF output for wind turbine model; minor fix in palmrungui

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