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

Last change on this file since 4275 was 4275, checked in by schwenkel, 4 years ago

move call of lpm at the end of intermediate timeloop and improve simple corrector particle interpolation method

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