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

Last change on this file since 3739 was 3739, checked in by dom_dwd_user, 3 years ago

biometeorology_mod.f90:
(N) Auto-adjusting thermal_comfort flag if not set by user, but thermal_indices set as output quantities.
(C) Renamed flags "bio_<index>" to "do_calculate_<index>" for better readability
(C) Removed everything related to "time_bio_results" as this is never used.
(C) Moved humidity warning to check_data_output so it will also be triggered if thermal_comofrt flag was auto-set later.
(B) Fixed bug in mrt calculation introduced with my commit yesterday.

time_integration.f90:
(C) Removed everything related to "time_bio_results" as this is never used

module_interface.f90:
(C) Removed call to empty method bio_check_parameters.

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