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

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