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

Last change on this file since 4180 was 4180, checked in by scharf, 5 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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