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

Last change on this file since 4180 was 4180, checked in by scharf, 2 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
Line 
1!> @file time_integration.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: time_integration.f90 4180 2019-08-21 14:37:54Z scharf $
27! copy diss, diss_p, tdiss_m to GPU
28!
29! 4144 2019-08-06 09:11:47Z raasch
30! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
31!
32! 4126 2019-07-30 11:09:11Z gronemeier
33! renamed routine to calculate uv exposure
34!
35! 4111 2019-07-22 18:16:57Z suehring
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
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
43! Moved call to radiation module out of intermediate time loop
44!
45! 4048 2019-06-21 21:00:21Z knoop
46! Moved production_e_init call into turbulence_closure_mod
47!
48! 4047 2019-06-21 18:58:09Z knoop
49! Added remainings of swap_timelevel upon its dissolution
50!
51! 4043 2019-06-18 16:59:00Z schwenkel
52! Further LPM modularization
53!
54! 4039 2019-06-18 10:32:41Z suehring
55! Rename subroutines in module for diagnostic quantities
56!
57! 4029 2019-06-14 14:04:35Z raasch
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
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
67! Mass (volume) flux correction included to ensure global mass conservation for child domains.
68!
69! 3994 2019-05-22 18:08:09Z suehring
70! output of turbulence intensity added
71!
72! 3988 2019-05-22 11:32:37Z kanani
73! Implement steerable output interval for virtual measurements
74!
75! 3968 2019-05-13 11:04:01Z suehring
76! replace nspec_out with n_matched_vars
77!
78! 3929 2019-04-24 12:52:08Z banzhafs
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
84! Changes related to global restructuring of location messages and introduction
85! of additional debug messages
86!
87! 3879 2019-04-08 20:25:23Z knoop
88! Moved wtm_forces to module_interface_actions
89!
90! 3872 2019-04-08 15:03:06Z knoop
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
101! added USE chem_gasphase_mod, replaced nspec by nspec since fixed compounds are not integrated
102!
103! 3820 2019-03-27 11:53:41Z forkel
104! renamed do_emiss to emissions_anthropogenic (ecc)
105!
106!
107! 3774 2019-03-04 10:52:49Z moh.hefny
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 $
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
116! Call indoor model after first timestep
117!
118! 3744 2019-02-15 18:38:58Z suehring
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
123! Removed everything related to "time_bio_results" as this is never used.
124!
125! 3724 2019-02-06 16:28:23Z kanani
126! Correct double-used log_point_s unit
127!
128! 3719 2019-02-06 13:10:18Z kanani
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
137! Data output for virtual measurements added
138!
139! 3704 2019-01-29 19:51:41Z suehring
140! Rename subroutines for surface-data output
141!
142! 3647 2019-01-02 14:10:44Z kanani
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!
147!
148! Description:
149! ------------
150!> Integration in time of the model equations, statistical analysis and graphic
151!> output
152!------------------------------------------------------------------------------!
153 SUBROUTINE time_integration
154 
155
156    USE advec_ws,                                                                                  &
157        ONLY:  ws_statistics
158
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
163
164    USE biometeorology_mod,                                                                        &
165        ONLY:  bio_calculate_thermal_index_maps, thermal_comfort, bio_calculate_uv_exposure,       &
166               uv_exposure
167
168    USE bulk_cloud_model_mod,                                                                      &
169        ONLY: bulk_cloud_model, calc_liquid_water_content, collision_turbulence,                   &
170              microphysics_morrison, microphysics_seifert
171
172    USE calc_mean_profile_mod,                                                                     &
173        ONLY:  calc_mean_profile
174
175    USE chem_emissions_mod,                                                                        &
176        ONLY:  chem_emissions_setup
177
178    USE chem_gasphase_mod,                                                                         &
179        ONLY:  nvar
180
181    USE chem_modules,                                                                              &
182        ONLY:  bc_cs_t_val, chem_species, cs_name, emissions_anthropogenic, n_matched_vars
183
184    USE chemistry_model_mod,                                                                       &
185        ONLY:  chem_boundary_conds
186
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,         &
196               dt_domask,dt_dopts, dt_dopr, dt_dopr_listing, dt_dots, dt_run_control,              &
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,    &
200               multi_agent_system_end, multi_agent_system_start, nesting_offline, neutral,         &
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,      &
209               time_do_sla, time_disturb, time_run_control, time_since_reference_point,            &
210               turbulent_inflow, turbulent_outflow, urban_surface,                                 &
211               use_initial_profile_as_reference, use_single_reference_value, u_gtrans, v_gtrans,   &
212               virtual_flight, virtual_measurement, ws_scheme_mom, ws_scheme_sca, timestep_count
213
214    USE cpulog,                                                                                    &
215        ONLY:  cpu_log, log_point, log_point_s
216
217    USE date_and_time_mod,                                                                         &
218        ONLY:  calc_date_and_time, hour_call_emis, hour_of_year
219
220    USE diagnostic_output_quantities_mod,                                                          &
221        ONLY:  doq_calculate,                                                                      &
222               timestep_number_at_prev_calc
223
224    USE flight_mod,                                                                                &
225        ONLY:  flight_measurement
226
227    USE indices,                                                                                   &
228        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, nzb, nzt
229
230    USE indoor_model_mod,                                                                          &
231        ONLY:  dt_indoor, im_main_heatcool, time_indoor
232
233    USE interfaces
234
235    USE kinds
236
237    USE land_surface_model_mod,                                                                    &
238        ONLY:  lsm_boundary_condition, lsm_energy_balance, lsm_soil_model, skip_time_do_lsm
239
240    USE lagrangian_particle_model_mod,                                                             &
241        ONLY:  lpm_data_output_ptseries, lpm_interaction_droplets_ptq
242
243    USE lsf_nudging_mod,                                                                           &
244        ONLY:  calc_tnudge, ls_forcing_surf, ls_forcing_vert, nudge_ref
245
246    USE module_interface,                                                                          &
247        ONLY:  module_interface_actions, module_interface_swap_timelevel
248
249    USE multi_agent_system_mod,                                                                    &
250        ONLY:  agents_active, multi_agent_system
251
252    USE nesting_offl_mod,                                                                          &
253        ONLY:  nesting_offl_bc, nesting_offl_mass_conservation
254
255    USE netcdf_data_input_mod,                                                                     &
256        ONLY:  chem_emis, chem_emis_att, nest_offl, netcdf_data_input_offline_nesting
257
258    USE ocean_mod,                                                                                 &
259        ONLY:  prho_reference
260
261    USE particle_attributes,                                                                       &
262        ONLY:  particle_advection, particle_advection_start, use_sgs_for_particles, wang_kernel
263
264    USE pegrid
265
266    USE pmc_interface,                                                                             &
267        ONLY:  nested_run, nesting_mode, pmci_boundary_conds, pmci_datatrans, pmci_synchronize,    &
268        pmci_ensure_nest_mass_conservation, pmci_ensure_nest_mass_conservation_vertical,           &
269        pmci_set_swaplevel
270
271    USE progress_bar,                                                                              &
272        ONLY:  finish_progress_bar, output_progress_bar
273
274    USE prognostic_equations_mod,                                                                  &
275        ONLY:  prognostic_equations_cache, prognostic_equations_vector
276
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
280
281    USE salsa_mod,                                                                                 &
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
285
286    USE spectra_mod,                                                                               &
287        ONLY: average_count_sp, averaging_interval_sp, calc_spectra, dt_dosp, skip_time_dosp
288
289    USE statistics,                                                                                &
290        ONLY:  flow_statistics_called, hom, pr_palm, sums_ls_l
291
292
293    USE surface_layer_fluxes_mod,                                                                  &
294        ONLY:  surface_layer_fluxes
295
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,               &
299               skip_time_dosurf_av, time_dosurf, time_dosurf_av
300
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,                                                                    &
309        ONLY:  tcm_diffusivities
310
311    USE urban_surface_mod,                                                                         &
312        ONLY:  usm_boundary_condition, usm_material_heat_model, usm_material_model,                &
313               usm_surface_energy_balance, usm_green_heat_model
314
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
319
320    USE virtual_measurement_mod,                                                                   &
321        ONLY:  dt_virtual_measurement,                                                             &
322               time_virtual_measurement,                                                           &
323               vm_data_output,                                                                     &
324               vm_sampling,                                                                        &
325               vm_time_start
326
327
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,        &
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
334
335    USE control_parameters,                                                                        &
336        ONLY:  tsc
337
338    USE indices,                                                                                   &
339        ONLY:  advc_flags_m, advc_flags_s, nyn, nyng, nys, nysg, nz, nzb_max, wall_flags_0
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, &
345               sums_wsnrs_ws_l, sums_wsss_ws_l, weight_substep, sums_salsa_ws_l
346
347    USE surface_mod,                                                                               &
348        ONLY:  bc_h, enter_surface_arrays, exit_surface_arrays
349#endif
350
351
352    IMPLICIT NONE
353
354    CHARACTER (LEN=9) ::  time_to_string   !<
355
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
364
365    REAL(wp) ::  dt_3d_old  !< temporary storage of timestep to be used for
366                            !< steering of run control output interval
367    REAL(wp) ::  time_since_reference_point_save  !< original value of
368                                                  !< time_since_reference_point
369
370
371!
372!-- Copy data from arrays_3d
373!$ACC DATA &
374!$ACC COPY(d(nzb+1:nzt,nys:nyn,nxl:nxr)) &
375!$ACC COPY(diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
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 &
386!$ACC COPY(diss_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
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)) &
393!$ACC COPY(tdiss_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
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)) &
407!$ACC COPYIN(heatflux_output_conversion(nzb:nzt+1)) &
408!$ACC COPYIN(momentumflux_output_conversion(nzb:nzt+1)) &
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
417!
418!-- Copy data from control_parameters
419!$ACC DATA &
420!$ACC COPYIN(tsc(1:5))
421
422!
423!-- Copy data from indices
424!$ACC DATA &
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)) &
427!$ACC COPYIN(wall_flags_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
428
429!
430!-- Copy data from surface_mod
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
440!
441!-- Copy data from statistics
442!$ACC DATA &
443!$ACC COPYIN(hom(0:nz+1,1:2,1:4,0)) &
444!$ACC COPYIN(rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions)) &
445!$ACC COPYIN(weight_substep(1:intermediate_timestep_count_max)) &
446!$ACC COPY(sums_l(nzb:nzt+1,1:pr_palm,0)) &
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
463#if defined( _OPENACC )
464    CALL enter_surface_arrays
465#endif
466
467!
468!-- At beginning determine the first time step
469    CALL timestep
470!
471!-- Synchronize the timestep in case of nested run.
472    IF ( nested_run )  THEN
473!
474!--    Synchronization by unifying the time step.
475!--    Global minimum of all time-steps is used for all.
476       CALL pmci_synchronize
477    ENDIF
478
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)
482!-- need to be determined first --> CALL flow_statistics at the beginning of
483!-- run_control
484    CALL run_control
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.
488    IF ( coupling_mode /= 'uncoupled'  .AND.  run_coupled  .AND. .NOT. vnested )  THEN
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.
492       DO WHILE ( time_coupling >= dt_coupling  .AND.  terminate_coupled == 0 )
493          CALL surface_coupler
494          time_coupling = time_coupling - dt_coupling
495       ENDDO
496       IF (time_coupling == 0.0_wp  .AND.  time_since_reference_point < dt_coupling )  THEN
497          time_coupling = time_since_reference_point
498       ENDIF
499    ENDIF
500
501    CALL location_message( 'atmosphere (and/or ocean) time-stepping', 'start' )
502
503!
504!-- Start of the time loop
505    DO  WHILE ( simulated_time < end_time  .AND.  .NOT. stop_dt  .AND. .NOT. terminate_run )
506
507       CALL cpu_log( log_point_s(10), 'timesteps', 'start' )
508!
509!--    Vertical nesting: initialize fine grid
510       IF ( vnested ) THEN
511          IF ( .NOT. vnest_init  .AND.  simulated_time >= vnest_start_time )  THEN
512             CALL cpu_log( log_point_s(22), 'vnest_init', 'start' )
513             CALL vnest_init_fine
514             vnest_init = .TRUE.
515             CALL cpu_log( log_point_s(22), 'vnest_init', 'stop' )
516          ENDIF
517       ENDIF
518!
519!--    Determine ug, vg and w_subs in dependence on data from external file
520!--    LSF_DATA
521       IF ( large_scale_forcing .AND. lsf_vert )  THEN
522           CALL ls_forcing_vert ( simulated_time )
523           sums_ls_l = 0.0_wp
524       ENDIF
525
526!
527!--    Set pt_init and q_init to the current profiles taken from
528!--    NUDGING_DATA
529       IF ( nudging )  THEN
530           CALL nudge_ref ( simulated_time )
531!
532!--        Store temperature gradient at the top boundary for possible Neumann
533!--        boundary condition
534           bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
535           bc_q_t_val  = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
536           IF ( air_chemistry )  THEN
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) )                            &
540                               / dzu(nzt+1)
541              ENDDO
542           ENDIF
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
560       ENDIF
561!
562!--    If forcing by larger-scale models is applied, check if new data
563!--    at domain boundaries need to be read.
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
567       ENDIF
568
569!
570!--    Execute all other module actions routunes
571       CALL module_interface_actions( 'before_timestep' )
572       
573!
574!--    Start of intermediate step loop
575       intermediate_timestep_count = 0
576       DO  WHILE ( intermediate_timestep_count < intermediate_timestep_count_max )
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!
586!--       Calculate those variables needed in the tendency terms which need
587!--       global communication
588          IF ( .NOT. use_single_reference_value  .AND.  .NOT. use_initial_profile_as_reference )   &
589          THEN
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!)
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
598             IF ( ocean_mode )  THEN
599                CALL calc_mean_profile( rho_ocean, 64 )
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
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 ).
610             ref_state = MERGE( MAXVAL(ref_state), ref_state, ref_state == 0.0_wp )
611          ENDIF
612
613          IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  intermediate_timestep_count == 1 )     &
614          THEN
615             CALL ws_statistics
616          ENDIF
617!
618!--       In case of nudging calculate current nudging time scale and horizontal
619!--       means of u, v, pt and q
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
627!
628!--       Execute all other module actions routunes
629          CALL module_interface_actions( 'before_prognostic_equations' )
630!
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.
636          IF ( loop_optimization == 'cache' )  THEN
637             CALL prognostic_equations_cache
638          ELSEIF ( loop_optimization == 'vector' )  THEN
639             CALL prognostic_equations_vector
640          ENDIF
641
642!
643!--       Particle transport/physics with the Lagrangian particle model
644!--       (only once during intermediate steps, because it uses an Euler-step)
645!--       ### particle model should be moved before prognostic_equations, in order
646!--       to regard droplet interactions directly
647          IF ( particle_advection  .AND.  time_since_reference_point >= particle_advection_start   &
648               .AND.  intermediate_timestep_count == 1 )                                           &
649          THEN
650             CALL module_interface_actions( 'after_prognostic_equations' )
651             first_call_lpm = .FALSE.
652          ENDIF
653
654!
655!--       Interaction of droplets with temperature and mixing ratio.
656!--       Droplet condensation and evaporation is calculated within
657!--       advec_particles.
658          IF ( cloud_droplets  .AND.  intermediate_timestep_count == intermediate_timestep_count_max ) &
659          THEN
660             CALL lpm_interaction_droplets_ptq
661          ENDIF
662
663!
664!--       Movement of agents in multi agent system
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
669             CALL multi_agent_system
670             first_call_mas = .FALSE.
671          ENDIF
672
673!
674!--       Exchange of ghost points (lateral boundary conditions)
675          CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'start' )
676
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 )
682          IF ( rans_tke_e  .OR.  wang_kernel  .OR.  collision_turbulence                           &
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
690          IF ( ocean_mode )  THEN
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 )
697             IF ( bulk_cloud_model .AND. microphysics_morrison )  THEN
698                CALL exchange_horiz( qc_p, nbgp )
699                CALL exchange_horiz( nc_p, nbgp )
700             ENDIF
701             IF ( bulk_cloud_model .AND. microphysics_seifert )  THEN
702                CALL exchange_horiz( qr_p, nbgp )
703                CALL exchange_horiz( nr_p, nbgp )
704             ENDIF
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
712          IF ( passive_scalar )  CALL exchange_horiz( s_p, nbgp )
713          IF ( air_chemistry )  THEN
714             DO  lsp = 1, nvar
715                CALL exchange_horiz( chem_species(lsp)%conc_p, nbgp )
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
735!
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
746                lsp_usr = 1
747                DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
748                   IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN
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
755
756          ENDIF
757
758          IF ( salsa  .AND.  time_since_reference_point >= skip_time_do_salsa )  THEN
759
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 )
765                ENDDO
766             ENDDO
767             IF ( .NOT. salsa_gases_from_chem )  THEN
768                DO  ig = 1, ngases_salsa
769                   CALL salsa_boundary_conds( salsa_gas(ig)%conc_p, salsa_gas(ig)%init )
770                ENDDO
771             ENDIF
772
773          ENDIF
774
775!
776!--       Incrementing timestep counter
777          timestep_count = timestep_count + 1
778
779          CALL cpu_log( log_point(28), 'swap_timelevel', 'start' )
780!
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!
790!--       Vertical nesting: Interpolate fine grid data to the coarse grid
791          IF ( vnest_init ) THEN
792             CALL cpu_log( log_point_s(37), 'vnest_anterpolate', 'start' )
793             CALL vnest_anterpolate
794             CALL cpu_log( log_point_s(37), 'vnest_anterpolate', 'stop' )
795          ENDIF
796
797          IF ( nested_run )  THEN
798
799             CALL cpu_log( log_point(60), 'nesting', 'start' )
800!
801!--          Domain nesting. The data transfer subroutines pmci_parent_datatrans
802!--          and pmci_child_datatrans are called inside the wrapper
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 )
807
808             IF ( TRIM( nesting_mode ) == 'two-way' .OR.  nesting_mode == 'vertical' )  THEN
809
810                CALL cpu_log( log_point_s(92), 'exchange-horiz-nest', 'start' )
811!
812!--             Exchange_horiz is needed for all parent-domains after the
813!--             anterpolation
814                CALL exchange_horiz( u, nbgp )
815                CALL exchange_horiz( v, nbgp )
816                CALL exchange_horiz( w, nbgp )
817                IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
818
819                IF ( humidity )  THEN
820
821                   CALL exchange_horiz( q, nbgp )
822
823                   IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
824                       CALL exchange_horiz( qc, nbgp )
825                       CALL exchange_horiz( nc, nbgp )
826                   ENDIF
827                   IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
828                       CALL exchange_horiz( qr, nbgp )
829                       CALL exchange_horiz( nr, nbgp )
830                   ENDIF
831
832                ENDIF
833
834                IF ( passive_scalar )  CALL exchange_horiz( s, nbgp ) 
835
836                IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e, nbgp )
837
838                IF ( .NOT. constant_diffusion  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
839                   CALL exchange_horiz( diss, nbgp )
840                ENDIF
841
842                IF ( air_chemistry )  THEN
843                   DO  n = 1, nvar
844                      CALL exchange_horiz( chem_species(n)%conc, nbgp ) 
845                   ENDDO
846                ENDIF
847
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
862                CALL cpu_log( log_point_s(92), 'exchange-horiz-nest', 'stop' )
863
864             ENDIF
865
866!
867!--          Set boundary conditions again after interpolation and anterpolation.
868             CALL pmci_boundary_conds
869
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
902             CALL cpu_log( log_point(60), 'nesting', 'stop' )
903
904          ENDIF
905
906!
907!--       Temperature offset must be imposed at cyclic boundaries in x-direction
908!--       when a sloping surface is used
909          IF ( sloping_surface )  THEN
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
912          ENDIF
913
914!
915!--       Impose a turbulent inflow using the recycling method
916          IF ( turbulent_inflow )  CALL inflow_turbulence
917
918!
919!--       Set values at outflow boundary using the special outflow condition
920          IF ( turbulent_outflow )  CALL outflow_turbulence
921
922!
923!--       Impose a random perturbation on the horizontal velocity field
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 ) ) &
927          THEN
928             time_disturb = time_disturb + dt_3d
929             IF ( time_disturb >= dt_disturb )  THEN
930                IF ( disturbance_energy_limit /= 0.0_wp  .AND.                                     &
931                     hom(nzb+5,1,pr_palm,0) < disturbance_energy_limit )  THEN
932                   CALL disturb_field( 'u', tend, u )
933                   CALL disturb_field( 'v', tend, v )
934                ELSEIF ( ( .NOT. bc_lr_cyc  .OR.  .NOT. bc_ns_cyc )                                &
935                         .AND. .NOT. child_domain  .AND.  .NOT.  nesting_offline )                 &
936                THEN
937!
938!--                Runs with a non-cyclic lateral wall need perturbations
939!--                near the inflow throughout the whole simulation
940                   dist_range = 1
941                   CALL disturb_field( 'u', tend, u )
942                   CALL disturb_field( 'v', tend, v )
943                   dist_range = 0
944                ENDIF
945                time_disturb = time_disturb - dt_disturb
946             ENDIF
947          ENDIF
948
949!
950!--       Map forcing data derived from larger scale model onto domain
951!--       boundaries.
952          IF ( nesting_offline  .AND.  intermediate_timestep_count ==                              &
953                                       intermediate_timestep_count_max  )                          &
954             CALL nesting_offl_bc
955!
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
959             CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
960             CALL stg_main
961             CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
962          ENDIF
963!
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!
970!--       Reduce the velocity divergence via the equation for perturbation
971!--       pressure.
972          IF ( intermediate_timestep_count == 1  .OR. &
973                call_psolver_at_all_substeps )  THEN
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
981                CALL cpu_log( log_point_s(30), 'vnest_bc', 'start' )
982                CALL vnest_boundary_conds
983                CALL cpu_log( log_point_s(30), 'vnest_bc', 'stop' )
984 
985                IF ( coupling_mode == 'vnested_fine' )  CALL pres
986
987!--             Anterpolate TKE, satisfy Germano Identity
988                CALL cpu_log( log_point_s(28), 'vnest_anter_e', 'start' )
989                CALL vnest_anterpolate_e
990                CALL cpu_log( log_point_s(28), 'vnest_anter_e', 'stop' )
991
992             ELSE
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               
1003                CALL pres
1004
1005             ENDIF
1006
1007          ENDIF
1008
1009!
1010!--       If required, compute liquid water content
1011          IF ( bulk_cloud_model )  THEN
1012             CALL calc_liquid_water_content
1013          ENDIF
1014!
1015!--       If required, compute virtual potential temperature
1016          IF ( humidity )  THEN
1017             CALL compute_vpt
1018          ENDIF
1019
1020!
1021!--       Compute the diffusion quantities
1022          IF ( .NOT. constant_diffusion )  THEN
1023
1024!
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
1028             IF ( ( large_scale_forcing .AND. lsf_surf ) .AND.                                     &
1029                 intermediate_timestep_count == intermediate_timestep_count_max )                  &
1030             THEN
1031                CALL ls_forcing_surf( simulated_time )
1032             ENDIF
1033
1034!
1035!--          First the vertical (and horizontal) fluxes in the surface
1036!--          (constant flux) layer are computed
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' )
1041             ENDIF
1042!
1043!--          If required, solve the energy balance for the surface and run soil
1044!--          model. Call for horizontal as well as vertical surfaces
1045             IF ( land_surface .AND. time_since_reference_point >= skip_time_do_lsm)  THEN
1046
1047                CALL cpu_log( log_point(54), 'land_surface', 'start' )
1048!
1049!--             Call for horizontal upward-facing surfaces
1050                CALL lsm_energy_balance( .TRUE., -1 )
1051                CALL lsm_soil_model( .TRUE., -1, .TRUE. )
1052!
1053!--             Call for northward-facing surfaces
1054                CALL lsm_energy_balance( .FALSE., 0 )
1055                CALL lsm_soil_model( .FALSE., 0, .TRUE. )
1056!
1057!--             Call for southward-facing surfaces
1058                CALL lsm_energy_balance( .FALSE., 1 )
1059                CALL lsm_soil_model( .FALSE., 1, .TRUE. )
1060!
1061!--             Call for eastward-facing surfaces
1062                CALL lsm_energy_balance( .FALSE., 2 )
1063                CALL lsm_soil_model( .FALSE., 2, .TRUE. )
1064!
1065!--             Call for westward-facing surfaces
1066                CALL lsm_energy_balance( .FALSE., 3 )
1067                CALL lsm_soil_model( .FALSE., 3, .TRUE. )
1068               
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
1074
1075               
1076                CALL cpu_log( log_point(54), 'land_surface', 'stop' )
1077             ENDIF
1078!
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' )
1083               
1084                CALL usm_surface_energy_balance( .FALSE. )
1085                IF ( usm_material_model )  THEN
1086                   CALL usm_green_heat_model
1087                   CALL usm_material_heat_model ( .FALSE. )
1088                ENDIF
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
1096                CALL cpu_log( log_point(74), 'urban_surface', 'stop' )
1097             ENDIF
1098!
1099!--          Compute the diffusion coefficients
1100             CALL cpu_log( log_point(17), 'diffusivities', 'start' )
1101             IF ( .NOT. humidity ) THEN
1102                IF ( ocean_mode )  THEN
1103                   CALL tcm_diffusivities( prho, prho_reference )
1104                ELSE
1105                   CALL tcm_diffusivities( pt, pt_reference )
1106                ENDIF
1107             ELSE
1108                CALL tcm_diffusivities( vpt, pt_reference )
1109             ENDIF
1110             CALL cpu_log( log_point(17), 'diffusivities', 'stop' )
1111!
1112!--          Vertical nesting: set fine grid eddy viscosity top boundary condition
1113             IF ( vnest_init )  CALL vnest_boundary_conds_khkm
1114
1115          ENDIF
1116
1117       ENDDO   ! Intermediate step loop
1118
1119!
1120!--    Will be used at some point by flow_statistics.
1121       !$ACC UPDATE &
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!
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!
1178!--    If required, consider chemical emissions
1179       IF ( air_chemistry  .AND.  emissions_anthropogenic )  THEN
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
1185          IF ( hour_of_year  >  hour_call_emis )  THEN
1186             CALL chem_emissions_setup( chem_emis_att, chem_emis, n_matched_vars )
1187             hour_call_emis = hour_of_year
1188          ENDIF
1189       ENDIF
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
1196
1197       ENDIF
1198!
1199!--    If required, calculate indoor temperature, waste heat, heat flux
1200!--    through wall, etc.
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.
1204       IF ( indoor_model )  THEN
1205
1206          time_indoor = time_indoor + dt_3d
1207
1208          IF ( time_indoor >= dt_indoor  .OR.  current_timestep_number == 0 )  THEN
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!
1219!--    Increase simulation time and output times
1220       nr_timesteps_this_run      = nr_timesteps_this_run + 1
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
1224       simulated_time_chr         = time_to_string( time_since_reference_point )
1225
1226
1227
1228
1229       IF ( time_since_reference_point >= skip_time_data_output_av )  THEN
1230          time_do_av         = time_do_av       + dt_3d
1231       ENDIF
1232       IF ( time_since_reference_point >= skip_time_do2d_xy )  THEN
1233          time_do2d_xy       = time_do2d_xy     + dt_3d
1234       ENDIF
1235       IF ( time_since_reference_point >= skip_time_do2d_xz )  THEN
1236          time_do2d_xz       = time_do2d_xz     + dt_3d
1237       ENDIF
1238       IF ( time_since_reference_point >= skip_time_do2d_yz )  THEN
1239          time_do2d_yz       = time_do2d_yz     + dt_3d
1240       ENDIF
1241       IF ( time_since_reference_point >= skip_time_do3d    )  THEN
1242          time_do3d          = time_do3d        + dt_3d
1243       ENDIF
1244       DO  mid = 1, masks
1245          IF ( time_since_reference_point >= skip_time_domask(mid) )  THEN
1246             time_domask(mid)= time_domask(mid) + dt_3d
1247          ENDIF
1248       ENDDO
1249       IF ( time_since_reference_point >= skip_time_dosp )  THEN
1250          time_dosp       = time_dosp        + dt_3d
1251       ENDIF
1252       time_dots          = time_dots        + dt_3d
1253       IF ( .NOT. first_call_lpm )  THEN
1254          time_dopts      = time_dopts       + dt_3d
1255       ENDIF
1256       IF ( time_since_reference_point >= skip_time_dopr )  THEN
1257          time_dopr       = time_dopr        + dt_3d
1258       ENDIF
1259       time_dopr_listing  = time_dopr_listing + dt_3d
1260       time_run_control   = time_run_control + dt_3d
1261!
1262!--    Increment time-counter for surface output
1263       IF ( surface_output )  THEN
1264          IF ( time_since_reference_point >= skip_time_dosurf )  THEN
1265             time_dosurf    = time_dosurf + dt_3d
1266          ENDIF
1267          IF ( time_since_reference_point >= skip_time_dosurf_av )  THEN
1268             time_dosurf_av = time_dosurf_av + dt_3d
1269          ENDIF
1270       ENDIF
1271!
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!
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
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
1288          ENDIF
1289          time_stg_call = time_stg_call + dt_3d
1290       ENDIF
1291
1292!
1293!--    Data exchange between coupled models
1294       IF ( coupling_mode /= 'uncoupled'  .AND.  run_coupled  .AND.  .NOT. vnested )  THEN
1295          time_coupling = time_coupling + dt_3d
1296
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.
1303          DO WHILE ( time_coupling >= dt_coupling  .AND.  terminate_coupled == 0 )
1304             CALL surface_coupler
1305             time_coupling = time_coupling - dt_coupling
1306          ENDDO
1307       ENDIF
1308
1309!
1310!--    Biometeorology calculation of stationary thermal indices
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!!
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 ) &
1319                    ) )  THEN
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
1328             CALL bio_calculate_uv_exposure
1329          ENDIF
1330       ENDIF
1331
1332!
1333!--    Execute alle other module actions routunes
1334       CALL module_interface_actions( 'after_integration' )
1335
1336!
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
1347!--    This call is skipped if the remote model has already initiated a restart.
1348       IF ( .NOT. terminate_run )  CALL check_for_restart
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
1363       IF ( averaging_interval_pr /= 0.0_wp  .AND.                                                 &
1364            ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.                                &
1365            time_since_reference_point >= skip_time_dopr )  THEN
1366          time_dopr_av = time_dopr_av + dt_3d
1367          IF ( time_dopr_av >= dt_averaging_input_pr )  THEN
1368             do_sum = .TRUE.
1369             time_dopr_av = MOD( time_dopr_av, MAX( dt_averaging_input_pr, dt_3d ) )
1370          ENDIF
1371       ENDIF
1372       IF ( do_sum )  CALL flow_statistics
1373
1374!
1375!--    Sum-up 3d-arrays for later output of time-averaged 2d/3d/masked data
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 )                               &
1379       THEN
1380          time_do_sla = time_do_sla + dt_3d
1381          IF ( time_do_sla >= dt_averaging_input )  THEN
1382             IF ( current_timestep_number > timestep_number_at_prev_calc )                         &
1383                CALL doq_calculate
1384
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
1390!
1391!--    Average surface data
1392       IF ( surface_output )  THEN
1393          IF ( averaging_interval_surf /= 0.0_wp                                                   &
1394                .AND.  ( dt_dosurf_av - time_dosurf_av ) <= averaging_interval_surf                &
1395                .AND.  time_since_reference_point >= skip_time_dosurf_av )  THEN
1396             IF ( time_dosurf_av >= dt_averaging_input )  THEN       
1397                CALL surface_data_output_averaging
1398                average_count_surf = average_count_surf + 1
1399             ENDIF
1400          ENDIF
1401       ENDIF
1402
1403!
1404!--    Calculate spectra for time averaging
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
1407          time_dosp_av = time_dosp_av + dt_3d
1408          IF ( time_dosp_av >= dt_averaging_input_pr )  THEN
1409             CALL calc_spectra
1410             time_dosp_av = MOD( time_dosp_av, MAX( dt_averaging_input_pr, dt_3d ) )
1411          ENDIF
1412       ENDIF
1413
1414!
1415!--    Call flight module and output data
1416       IF ( virtual_flight )  THEN
1417          CALL flight_measurement
1418          CALL data_output_flight
1419       ENDIF
1420!
1421!--    Take virtual measurements
1422       IF ( virtual_measurement  .AND.  time_virtual_measurement >= dt_virtual_measurement         &
1423                                 .AND.  vm_time_start <= time_since_reference_point )  THEN
1424          CALL vm_sampling
1425          CALL vm_data_output
1426          time_virtual_measurement = MOD(      time_virtual_measurement,                           &
1427                                          MAX( dt_virtual_measurement, dt_3d ) )
1428       ENDIF
1429!
1430!--    Profile output (ASCII) on file
1431       IF ( time_dopr_listing >= dt_dopr_listing )  THEN
1432          CALL print_1d
1433          time_dopr_listing = MOD( time_dopr_listing, MAX( dt_dopr_listing, dt_3d ) )
1434       ENDIF
1435
1436!
1437!--    Graphic output for PROFIL
1438       IF ( time_dopr >= dt_dopr  .AND.  time_since_reference_point >= skip_time_dopr )  THEN
1439          IF ( dopr_n /= 0 )  CALL data_output_profiles
1440          time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
1441          time_dopr_av = 0.0_wp    ! due to averaging (see above)
1442       ENDIF
1443
1444!
1445!--    Graphic output for time series
1446       IF ( time_dots >= dt_dots )  THEN
1447          CALL data_output_tseries
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
1454       IF ( time_dosp >= dt_dosp  .AND.  time_since_reference_point >= skip_time_dosp )  THEN
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)
1462       IF ( time_do2d_xy >= dt_do2d_xy  .AND.  time_since_reference_point >= skip_time_do2d_xy )  THEN
1463          IF ( current_timestep_number > timestep_number_at_prev_calc )                            &
1464             CALL doq_calculate
1465
1466          CALL data_output_2d( 'xy', 0 )
1467          time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
1468       ENDIF
1469       IF ( time_do2d_xz >= dt_do2d_xz  .AND.  time_since_reference_point >= skip_time_do2d_xz )  THEN
1470          IF ( current_timestep_number > timestep_number_at_prev_calc )                            &
1471
1472             CALL doq_calculate
1473          CALL data_output_2d( 'xz', 0 )
1474          time_do2d_xz = MOD( time_do2d_xz, MAX( dt_do2d_xz, dt_3d ) )
1475       ENDIF
1476       IF ( time_do2d_yz >= dt_do2d_yz  .AND.  time_since_reference_point >= skip_time_do2d_yz )  THEN
1477          IF ( current_timestep_number > timestep_number_at_prev_calc )                            &
1478             CALL doq_calculate
1479
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)
1486       IF ( time_do3d >= dt_do3d  .AND.  time_since_reference_point >= skip_time_do3d )  THEN
1487          IF ( current_timestep_number > timestep_number_at_prev_calc )                            &
1488             CALL doq_calculate
1489
1490          CALL data_output_3d( 0 )
1491          time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
1492       ENDIF
1493
1494!
1495!--    Masked data output
1496       DO  mid = 1, masks
1497          IF ( time_domask(mid) >= dt_domask(mid)                                                  &
1498               .AND.  time_since_reference_point >= skip_time_domask(mid) )  THEN
1499             IF ( current_timestep_number > timestep_number_at_prev_calc )                         &
1500                CALL doq_calculate
1501
1502             CALL data_output_mask( 0, mid )
1503             time_domask(mid) = MOD( time_domask(mid), MAX( dt_domask(mid), dt_3d ) )
1504          ENDIF
1505       ENDDO
1506
1507!
1508!--    Output of time-averaged 2d/3d/masked data
1509       IF ( time_do_av >= dt_data_output_av                                                        &
1510            .AND.  time_since_reference_point >= skip_time_data_output_av )  THEN
1511          CALL average_3d_data
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
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 )
1521          DO  mid = 1, masks
1522             CALL data_output_mask( 1, mid )
1523          ENDDO
1524          time_do_av = MOD( time_do_av, MAX( dt_data_output_av, dt_3d ) )
1525       ENDIF
1526!
1527!--    Output of surface data, instantaneous and averaged data
1528       IF ( surface_output )  THEN
1529          IF ( time_dosurf >= dt_dosurf  .AND.  time_since_reference_point >= skip_time_dosurf )  THEN
1530             CALL surface_data_output( 0 )
1531             time_dosurf = MOD( time_dosurf, MAX( dt_dosurf, dt_3d ) )
1532          ENDIF
1533          IF ( time_dosurf_av >= dt_dosurf_av  .AND.  time_since_reference_point >= skip_time_dosurf_av )  THEN
1534             CALL surface_data_output( 1 )
1535             time_dosurf_av = MOD( time_dosurf_av, MAX( dt_dosurf_av, dt_3d ) )
1536          ENDIF
1537       ENDIF
1538
1539!
1540!--    Output of particle time series
1541       IF ( particle_advection )  THEN
1542          IF ( time_dopts >= dt_dopts  .OR.                                                        &
1543               ( time_since_reference_point >= particle_advection_start  .AND.                     &
1544                 first_call_lpm ) )  THEN
1545             CALL lpm_data_output_ptseries
1546             time_dopts = MOD( time_dopts, MAX( dt_dopts, dt_3d ) )
1547          ENDIF
1548       ENDIF
1549
1550!
1551!--    If required, set the heat flux for the next time step to a random value
1552       IF ( constant_heatflux  .AND.  random_heatflux )  THEN
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
1568       ENDIF
1569
1570!
1571!--    Execute alle other module actions routunes
1572       CALL module_interface_actions( 'after_timestep' )
1573
1574!
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!
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!
1591!--    Computation and output of run control parameters.
1592!--    This is also done whenever perturbations have been imposed
1593       IF ( time_run_control >= dt_run_control  .OR.                                               &
1594            timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created )                           &
1595       THEN
1596          CALL run_control
1597          IF ( time_run_control >= dt_run_control )  THEN
1598             time_run_control = MOD( time_run_control, MAX( dt_run_control, dt_3d_old ) )
1599          ENDIF
1600       ENDIF
1601
1602!
1603!--    Output elapsed simulated time in form of a progress bar on stdout
1604       IF ( myid == 0 )  CALL output_progress_bar
1605
1606       CALL cpu_log( log_point_s(10), 'timesteps', 'stop' )
1607
1608
1609    ENDDO   ! time loop
1610
1611#if defined( _OPENACC )
1612    CALL exit_surface_arrays
1613#endif
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
1622!
1623!-- Vertical nesting: Deallocate variables initialized for vertical nesting   
1624    IF ( vnest_init )  CALL vnest_deallocate
1625
1626    IF ( myid == 0 )  CALL finish_progress_bar
1627
1628    CALL location_message( 'atmosphere (and/or ocean) time-stepping', 'finished' )
1629
1630 END SUBROUTINE time_integration
Note: See TracBrowser for help on using the repository browser.