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

Last change on this file since 2865 was 2817, checked in by knoop, 7 years ago

Preliminary gust module interface implemented

  • Property svn:keywords set to Id
File size: 50.9 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!
[2718]17! Copyright 1997-2018 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[484]20! Current revisions:
[1092]21! ------------------
[1928]22!
[2233]23!
[1366]24! Former revisions:
25! -----------------
26! $Id: time_integration.f90 2817 2018-02-19 16:32:21Z kanani $
[2817]27! Preliminary gust module interface implemented
28!
29! 2801 2018-02-14 16:01:55Z thiele
[2801]30! Changed lpm from subroutine to module.
31! Introduce particle transfer in nested models.
32!
33! 2776 2018-01-31 10:44:42Z Giersch
[2776]34! Variable use_synthetic_turbulence_generator has been abbreviated
35!
36! 2773 2018-01-30 14:12:54Z suehring
[2773]37! - Nesting for chemical species
38!
39! 2766 2018-01-22 17:17:47Z kanani
[2766]40! Removed preprocessor directive __chem
41!
42! 2718 2018-01-02 08:49:38Z maronga
[2716]43! Corrected "Former revisions" section
44!
45! 2696 2017-12-14 17:12:51Z kanani
46! - Change in file header (GPL part)
[2696]47! - Implementation of uv exposure model (FK)
48! - Moved vnest_boundary_conds_khkm from tcm_diffusivities to here (TG)
49! - renamed diffusivities to tcm_diffusivities (TG)
50! - implement prognostic equation for diss (TG)
51! - Moved/commented CALL to chem_emissions (FK)
52! - Added CALL to chem_emissions (FK)
53! - Implementation of chemistry module (FK)
54! - Calls for setting boundary conditions in USM and LSM (MS)
55! - Large-scale forcing with larger-scale models implemented (MS)
56! - Rename usm_radiation into radiation_interactions; merge with branch
57!   radiation (MS)
58! - added call for usm_green_heat_model for green building surfaces (RvT)
59! - added call for usm_temperature_near_surface for use in indoor model (RvT)
60!
61! 2617 2017-11-16 12:47:24Z suehring
[2617]62! Bugfix, assure that the reference state does not become zero.
63!
64! 2563 2017-10-19 15:36:10Z Giersch
[2563]65! Variable wind_turbine moved to module control_parameters
66!
67! 2365 2017-08-21 14:59:59Z kanani
[2365]68! Vertical grid nesting implemented (SadiqHuq)
69!
70! 2320 2017-07-21 12:47:43Z suehring
[2311]71! Set bottom boundary conditions after nesting interpolation and anterpolation
72!
73! 2299 2017-06-29 10:14:38Z maronga
[2299]74! Call of soil model adjusted
75!
76! 2292 2017-06-20 09:51:42Z schwenkel
[2292]77! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
78! includes two more prognostic equations for cloud drop concentration (nc) 
79! and cloud water content (qc).
80!
81! 2271 2017-06-09 12:34:55Z sward
[2271]82! Start timestep message changed
83!
84! 2259 2017-06-08 09:09:11Z gronemeier
[2259]85! Implemented synthetic turbulence generator
[1366]86!
[2259]87! 2233 2017-05-30 18:08:54Z suehring
88!
[2233]89! 2232 2017-05-30 17:47:52Z suehring
90! Adjustments to new topography and surface concept
91! Modify passed parameters for disturb_field
92!
[2179]93! 2178 2017-03-17 11:07:39Z hellstea
94! Setting perturbations at all times near inflow boundary is removed
95! in case of nested boundaries
96!
[2175]97! 2174 2017-03-13 08:18:57Z maronga
98! Added support for nesting with cloud microphysics
99!
[2119]100! 2118 2017-01-17 16:38:49Z raasch
101! OpenACC directives and related code removed
102!
[2051]103! 2050 2016-11-08 15:00:55Z gronemeier
104! Implement turbulent outflow condition
105!
[2032]106! 2031 2016-10-21 15:11:58Z knoop
107! renamed variable rho to rho_ocean
108!
[2012]109! 2011 2016-09-19 17:29:57Z kanani
110! Flag urban_surface is now defined in module control_parameters,
111! removed commented CALLs of global_min_max.
112!
[2008]113! 2007 2016-08-24 15:47:17Z kanani
114! Added CALLs for new urban surface model
115!
[2001]116! 2000 2016-08-20 18:09:15Z knoop
117! Forced header and separation lines into 80 columns
118!
[1977]119! 1976 2016-07-27 13:28:04Z maronga
120! Simplified calls to radiation model
121!
[1961]122! 1960 2016-07-12 16:34:24Z suehring
123! Separate humidity and passive scalar
124!
[1958]125! 1957 2016-07-07 10:43:48Z suehring
126! flight module added
127!
[1933]128! 1919 2016-05-27 14:51:23Z raasch
129! Initial version of purely vertical nesting introduced.
[1928]130!
[1919]131! 1918 2016-05-27 14:35:57Z raasch
132! determination of time step moved to the end of the time step loop,
133! the first time step is now always calculated before the time step loop (i.e.
134! also in case of restart runs)
135!
[1917]136! 1914 2016-05-26 14:44:07Z witha
137! Added call for wind turbine model
138!
[1879]139! 1878 2016-04-19 12:30:36Z hellstea
140! Synchronization for nested runs rewritten
141!
[1854]142! 1853 2016-04-11 09:00:35Z maronga
143! Adjusted for use with radiation_scheme = constant
144!
[1851]145! 1849 2016-04-08 11:33:18Z hoffmann
146! Adapted for modularization of microphysics
147!
[1834]148! 1833 2016-04-07 14:23:03Z raasch
149! spectrum renamed spectra_mod, spectra related variables moved to spectra_mod
150!
[1832]151! 1831 2016-04-07 13:15:51Z hoffmann
152! turbulence renamed collision_turbulence
153!
[1823]154! 1822 2016-04-07 07:49:42Z hoffmann
155! icloud_scheme replaced by microphysics_*
156!
[1809]157! 1808 2016-04-05 19:44:00Z raasch
158! output message in case unscheduled radiation calls removed
159!
[1798]160! 1797 2016-03-21 16:50:28Z raasch
161! introduction of different datatransfer modes
162!
[1792]163! 1791 2016-03-11 10:41:25Z raasch
164! call of pmci_update_new removed
165!
[1787]166! 1786 2016-03-08 05:49:27Z raasch
167! +module spectrum
168!
[1784]169! 1783 2016-03-06 18:36:17Z raasch
170! switch back of netcdf data format for mask output moved to the mask output
171! routine
172!
[1782]173! 1781 2016-03-03 15:12:23Z raasch
174! some pmc calls removed at the beginning (before timeloop),
175! pmc initialization moved to the main program
176!
[1765]177! 1764 2016-02-28 12:45:19Z raasch
178! PMC_ACTIVE flags removed,
179! bugfix: nest synchronization after first call of timestep
180!
[1763]181! 1762 2016-02-25 12:31:13Z hellstea
182! Introduction of nested domain feature
183!
[1737]184! 1736 2015-12-04 08:56:33Z raasch
185! no perturbations added to total domain if energy limit has been set zero
186!
[1692]187! 1691 2015-10-26 16:17:44Z maronga
188! Added option for spin-ups without land surface and radiation models. Moved calls
189! for radiation and lan surface schemes.
190!
[1683]191! 1682 2015-10-07 23:56:08Z knoop
192! Code annotations made doxygen readable
193!
[1672]194! 1671 2015-09-25 03:29:37Z raasch
195! bugfix: ghostpoint exchange for array diss in case that sgs velocities are used
196! for particles
197!
[1586]198! 1585 2015-04-30 07:05:52Z maronga
199! Moved call of radiation scheme. Added support for RRTM
200!
[1552]201! 1551 2015-03-03 14:18:16Z maronga
202! Added interface for different radiation schemes.
203!
[1497]204! 1496 2014-12-02 17:25:50Z maronga
205! Added calls for the land surface model and radiation scheme
206!
[1403]207! 1402 2014-05-09 14:25:13Z raasch
208! location messages modified
209!
[1385]210! 1384 2014-05-02 14:31:06Z raasch
211! location messages added
212!
[1381]213! 1380 2014-04-28 12:40:45Z heinze
214! CALL of nudge_ref added
215! bc_pt_t_val and bc_q_t_val are updated in case nudging is used
216!
[1366]217! 1365 2014-04-22 15:03:56Z boeske
[1365]218! Reset sums_ls_l to zero at each timestep
219! +sums_ls_l
220! Calculation of reference state (previously in subroutine calc_mean_profile)
221
[1343]222! 1342 2014-03-26 17:04:47Z kanani
223! REAL constants defined as wp-kind
224!
[1321]225! 1320 2014-03-20 08:40:49Z raasch
[1320]226! ONLY-attribute added to USE-statements,
227! kind-parameters added to all INTEGER and REAL declaration statements,
228! kinds are defined in new module kinds,
229! old module precision_kind is removed,
230! revision history before 2012 removed,
231! comment fields (!:) to be used for variable explanations added to
232! all variable declaration statements
[1319]233! 1318 2014-03-17 13:35:16Z raasch
234! module interfaces removed
235!
[1309]236! 1308 2014-03-13 14:58:42Z fricke
237! +netcdf_data_format_save
238! For masked data, parallel netcdf output is not tested so far, hence
239! netcdf_data_format is switched back to non-paralell output.
240!
[1277]241! 1276 2014-01-15 13:40:41Z heinze
242! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
243!
[1258]244! 1257 2013-11-08 15:18:40Z raasch
245! acc-update-host directive for timestep removed
246!
[1242]247! 1241 2013-10-30 11:36:58Z heinze
248! Generalize calc_mean_profile for wider use
249! Determine shf and qsws in dependence on data from LSF_DATA
250! Determine ug and vg in dependence on data from LSF_DATA
[1222]251! 1221 2013-09-10 08:59:13Z raasch
252! host update of arrays before timestep is called
253!
[1182]254! 1179 2013-06-14 05:57:58Z raasch
255! mean profiles for reference state are only calculated if required,
256! small bugfix for background communication
257!
[1172]258! 1171 2013-05-30 11:27:45Z raasch
259! split of prognostic_equations deactivated (comment lines), for the time being
260!
[1132]261! 1128 2013-04-12 06:19:32Z raasch
[1128]262! asynchronous transfer of ghost point data realized for acc-optimized version:
263! prognostic_equations are first called two times for those points required for
264! the left-right and north-south exchange, respectively, and then for the
265! remaining points,
266! those parts requiring global communication moved from prognostic_equations to
267! here
[392]268!
[1116]269! 1115 2013-03-26 18:16:16Z hoffmann
270! calculation of qr and nr is restricted to precipitation
271!
[1114]272! 1113 2013-03-10 02:48:14Z raasch
273! GPU-porting of boundary conditions,
274! openACC directives updated
275! formal parameter removed from routine boundary_conds
276!
[1112]277! 1111 2013-03-08 23:54:10Z raasch
278! +internal timestep counter for cpu statistics added,
279! openACC directives updated
280!
[1093]281! 1092 2013-02-02 11:24:22Z raasch
282! unused variables removed
283!
[1066]284! 1065 2012-11-22 17:42:36Z hoffmann
285! exchange of diss (dissipation rate) in case of turbulence = .TRUE. added
286!
[1054]287! 1053 2012-11-13 17:11:03Z hoffmann
288! exchange of ghost points for nr, qr added
289!
[1037]290! 1036 2012-10-22 13:43:42Z raasch
291! code put under GPL (PALM 3.9)
292!
[1020]293! 1019 2012-09-28 06:46:45Z raasch
294! non-optimized version of prognostic_equations removed
295!
[1017]296! 1015 2012-09-27 09:23:24Z raasch
297! +call of prognostic_equations_acc
298!
[1002]299! 1001 2012-09-13 14:08:46Z raasch
300! all actions concerning leapfrog- and upstream-spline-scheme removed
301!
[850]302! 849 2012-03-15 10:35:09Z raasch
303! advec_particles renamed lpm, first_call_advec_particles renamed first_call_lpm
304!
[826]305! 825 2012-02-19 03:03:44Z raasch
306! wang_collision_kernel renamed wang_kernel
307!
[1]308! Revision 1.1  1997/08/11 06:19:04  raasch
309! Initial revision
310!
311!
312! Description:
313! ------------
[1682]314!> Integration in time of the model equations, statistical analysis and graphic
315!> output
[1]316!------------------------------------------------------------------------------!
[1682]317 SUBROUTINE time_integration
318 
[1]319
[1320]320    USE advec_ws,                                                              &
321        ONLY:  ws_statistics
322
323    USE arrays_3d,                                                             &
[2696]324        ONLY:  diss, diss_p, dzu, e, e_p, nc, nc_p, nr, nr_p, prho, pt, pt_p, pt_init, &
[2292]325               q_init, q, qc, qc_p, ql, ql_c, ql_v, ql_vp, qr, qr_p, q_p,      &
326               ref_state, rho_ocean, s, s_p, sa_p, tend, u, u_p, v, vpt,       &
327               v_p, w, w_p
[1320]328
[1365]329    USE calc_mean_profile_mod,                                                 &
[1320]330        ONLY:  calc_mean_profile
331
[2696]332    USE chemistry_model_mod,                                                   &
333        ONLY:  chem_emissions, chem_species
334
335    USE chem_modules,                                                          &
336        ONLY:  nspec 
337
[1320]338    USE control_parameters,                                                    &
[2696]339        ONLY:  advected_distance_x, advected_distance_y, air_chemistry,        &
340               average_count_3d, averaging_interval, averaging_interval_pr,    &
[1833]341               bc_lr_cyc, bc_ns_cyc, bc_pt_t_val,                              &
[1380]342               bc_q_t_val, call_psolver_at_all_substeps, cloud_droplets,       &
[1691]343               cloud_physics, constant_flux_layer, constant_heatflux,          &
344               create_disturbances, dopr_n, constant_diffusion, coupling_mode, &
345               coupling_start_time, current_timestep_number,                   &
346               disturbance_created, disturbance_energy_limit, dist_range,      &
347               do_sum, dt_3d, dt_averaging_input, dt_averaging_input_pr,       &
348               dt_coupling, dt_data_output_av, dt_disturb, dt_do2d_xy,         &
349               dt_do2d_xz, dt_do2d_yz, dt_do3d, dt_domask,dt_dopts, dt_dopr,   &
[2696]350               dt_dopr_listing, dt_dots, dt_dvrp, dt_run_control, end_time,    &
351               first_call_lpm, forcing, galilei_transformation, humidity,      &
[2232]352               intermediate_timestep_count, intermediate_timestep_count_max,   &
353               land_surface, large_scale_forcing,                              &
[1822]354               loop_optimization, lsf_surf, lsf_vert, masks,                   &
[2292]355               microphysics_morrison, microphysics_seifert, mid, nest_domain,  &
[1783]356               neutral, nr_timesteps_this_run, nudging,                        &
[2118]357               ocean, passive_scalar,                                          &
[1320]358               prho_reference, pt_reference, pt_slope_offset, random_heatflux, &
[2696]359               rans_tke_e, run_coupled, simulated_time, simulated_time_chr,    &
[1320]360               skip_time_do2d_xy, skip_time_do2d_xz, skip_time_do2d_yz,        &
361               skip_time_do3d, skip_time_domask, skip_time_dopr,               &
[1833]362               skip_time_data_output_av, sloping_surface,                      &
[1320]363               stop_dt, terminate_coupled, terminate_run, timestep_scheme,     &
364               time_coupling, time_do2d_xy, time_do2d_xz, time_do2d_yz,        &
365               time_do3d, time_domask, time_dopr, time_dopr_av,                &
366               time_dopr_listing, time_dopts, time_dosp, time_dosp_av,         &
367               time_dots, time_do_av, time_do_sla, time_disturb, time_dvrp,    &
[1496]368               time_run_control, time_since_reference_point,                   &
[2050]369               turbulent_inflow, turbulent_outflow, urban_surface,             &
[2011]370               use_initial_profile_as_reference,                               &
[2696]371               use_single_reference_value, uv_exposure, u_gtrans, v_gtrans,    &
372               virtual_flight, wind_turbine, ws_scheme_mom, ws_scheme_sca
[1320]373
374    USE cpulog,                                                                &
375        ONLY:  cpu_log, log_point, log_point_s
376
[1957]377    USE flight_mod,                                                            &
378        ONLY:  flight_measurement
379
[2817]380    USE gust_mod,                                                              &
381        ONLY:  gust_actions, gust_module_enabled
[1957]382
[1320]383    USE indices,                                                               &
[2232]384        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
[1320]385
386    USE interaction_droplets_ptq_mod,                                          &
387        ONLY:  interaction_droplets_ptq
388
[1918]389    USE interfaces
390
[1320]391    USE kinds
392
[1496]393    USE land_surface_model_mod,                                                &
[2696]394        ONLY:  lsm_boundary_condition, lsm_energy_balance, lsm_soil_model,     &
[1691]395               skip_time_do_lsm
[1496]396
[2320]397    USE lsf_nudging_mod,                                                       &
[2696]398        ONLY:  calc_tnudge, ls_forcing_surf, ls_forcing_vert, nudge_ref,       &
399               forcing_bc, forcing_bc_mass_conservation
[1320]400
[2696]401    USE netcdf_data_input_mod,                                                 &
402        ONLY:  force, netcdf_data_input_lsf
403
[1849]404    USE microphysics_mod,                                                      &
405        ONLY: collision_turbulence
406
[1320]407    USE particle_attributes,                                                   &
[1671]408        ONLY:  particle_advection, particle_advection_start,                   &
409               use_sgs_for_particles, wang_kernel
[1320]410
[1]411    USE pegrid
412
[1762]413    USE pmc_interface,                                                         &
[2311]414        ONLY:  nested_run, nesting_mode, pmci_boundary_conds, pmci_datatrans,  &
[1933]415               pmci_ensure_nest_mass_conservation, pmci_synchronize
[1762]416
[1402]417    USE progress_bar,                                                          &
418        ONLY:  finish_progress_bar, output_progress_bar
419
[1320]420    USE prognostic_equations_mod,                                              &
[2118]421        ONLY:  prognostic_equations_cache, prognostic_equations_vector
[1320]422
[1496]423    USE radiation_model_mod,                                                   &
[1976]424        ONLY: dt_radiation, force_radiation_call, radiation, radiation_control,&
[2696]425              radiation_interaction, radiation_interactions,                   &
[1976]426              skip_time_do_radiation, time_radiation
[1496]427
[1833]428    USE spectra_mod,                                                           &
429        ONLY: average_count_sp, averaging_interval_sp, calc_spectra, dt_dosp,  &
430              skip_time_dosp
[1786]431
[1320]432    USE statistics,                                                            &
[1918]433        ONLY:  flow_statistics_called, hom, pr_palm, sums_ls_l, u_max,         &
434               u_max_ijk, v_max, v_max_ijk, w_max, w_max_ijk
[1320]435
[1691]436    USE surface_layer_fluxes_mod,                                              &
437        ONLY:  surface_layer_fluxes
438
[2232]439    USE surface_mod,                                                           &
440        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
441
[2696]442    USE turbulence_closure_mod,                                                &
443        ONLY:  tcm_diffusivities, production_e_init
444
[2007]445    USE urban_surface_mod,                                                     &
[2696]446        ONLY:  usm_boundary_condition, usm_material_heat_model,                &
447               usm_material_model,                                             &
448               usm_surface_energy_balance, usm_green_heat_model,               &
449               usm_temperature_near_surface
[2007]450
[2259]451    USE synthetic_turbulence_generator_mod,                                    &
[2776]452        ONLY:  stg_main, use_syn_turb_gen
[2259]453
[1320]454    USE user_actions_mod,                                                      &
455        ONLY:  user_actions
456
[2696]457    USE uv_exposure_model_mod,                                                 &
458        ONLY:  uvem_calc_exposure
459
[1914]460    USE wind_turbine_model_mod,                                                &
[2563]461        ONLY:  wtm_forces
[1914]462
[2801]463    USE lpm_mod,                                                               &
464        ONLY:  lpm
465
[2365]466    USE vertical_nesting_mod,                                                  &
467        ONLY:  vnested, vnest_anterpolate, vnest_anterpolate_e,                &
468               vnest_boundary_conds, vnest_boundary_conds_khkm,                & 
469               vnest_deallocate, vnest_init, vnest_init_fine,                  &
470               vnest_start_time
471
[1]472    IMPLICIT NONE
473
[1682]474    CHARACTER (LEN=9) ::  time_to_string          !<
[2696]475    INTEGER           ::  n
476    INTEGER           ::  lsp
[1]477
[1918]478    REAL(wp) ::  dt_3d_old  !< temporary storage of timestep to be used for
479                            !< steering of run control output interval
[1]480!
[1918]481!-- At beginning determine the first time step
482    CALL timestep
[667]483
[1764]484!
485!-- Synchronize the timestep in case of nested run.
486    IF ( nested_run )  THEN
[1878]487!
488!--    Synchronization by unifying the time step.
489!--    Global minimum of all time-steps is used for all.
490       CALL pmci_synchronize
[1764]491    ENDIF
492
[1918]493!
494!-- Determine and print out the run control quantities before the first time
495!-- step of this run. For the initial run, some statistics (e.g. divergence)
496!-- need to be determined first.
497    IF ( simulated_time == 0.0_wp )  CALL flow_statistics
[1]498    CALL run_control
499
[108]500!
501!-- Data exchange between coupled models in case that a call has been omitted
502!-- at the end of the previous run of a job chain.
[2365]503    IF ( coupling_mode /= 'uncoupled'  .AND.  run_coupled .AND. .NOT. vnested)  THEN
[108]504!
505!--    In case of model termination initiated by the local model the coupler
506!--    must not be called because this would again cause an MPI hang.
[1918]507       DO WHILE ( time_coupling >= dt_coupling  .AND.  terminate_coupled == 0 )
[108]508          CALL surface_coupler
509          time_coupling = time_coupling - dt_coupling
510       ENDDO
[1918]511       IF (time_coupling == 0.0_wp  .AND.                                      &
512           time_since_reference_point < dt_coupling )                          &
[348]513       THEN
514          time_coupling = time_since_reference_point
515       ENDIF
[108]516    ENDIF
517
[1]518#if defined( __dvrp_graphics )
519!
520!-- Time measurement with dvrp software 
521    CALL DVRP_LOG_EVENT( 2, current_timestep_number )
522#endif
523
[2271]524    CALL location_message( 'starting timestep-sequence', .TRUE. )
[1]525!
526!-- Start of the time loop
527    DO  WHILE ( simulated_time < end_time  .AND.  .NOT. stop_dt  .AND. &
528                .NOT. terminate_run )
529
530       CALL cpu_log( log_point_s(10), 'timesteps', 'start' )
531!
[2365]532!--    Vertical nesting: initialize fine grid
533       IF ( vnested ) THEN
534          IF ( .NOT. vnest_init  .AND.  simulated_time >= vnest_start_time )  THEN
535             CALL cpu_log( log_point(80), 'vnest_init', 'start' )
536             CALL vnest_init_fine
537             vnest_init = .TRUE.
538             CALL cpu_log( log_point(80), 'vnest_init', 'stop' )
539          ENDIF
540       ENDIF
541!
[1241]542!--    Determine ug, vg and w_subs in dependence on data from external file
543!--    LSF_DATA
[1365]544       IF ( large_scale_forcing .AND. lsf_vert )  THEN
[1241]545           CALL ls_forcing_vert ( simulated_time )
[1365]546           sums_ls_l = 0.0_wp
[1241]547       ENDIF
548
549!
[1380]550!--    Set pt_init and q_init to the current profiles taken from
551!--    NUDGING_DATA
552       IF ( nudging )  THEN
553           CALL nudge_ref ( simulated_time )
554!
555!--        Store temperature gradient at the top boundary for possible Neumann
556!--        boundary condition
557           bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
558           bc_q_t_val  = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
559       ENDIF
[2696]560!
561!--    If forcing by larger-scale models is applied, check if new data
562!--    at domain boundaries need to be read.
563       IF ( forcing )  THEN
564          IF ( force%time(force%tind_p) <= simulated_time )                    &
565             CALL netcdf_data_input_lsf
566       ENDIF
[1380]567
568!
[2817]569!--    Execute the gust module actions
570       IF ( gust_module_enabled )  THEN
571          CALL gust_actions( 'before_timestep' )
572       ENDIF
573
574!
[1]575!--    Execute the user-defined actions
576       CALL user_actions( 'before_timestep' )
577
578!
[1914]579!--    Calculate forces by wind turbines
580       IF ( wind_turbine )  THEN
581
582          CALL cpu_log( log_point(55), 'wind_turbine', 'start' )
583
584          CALL wtm_forces
585
586          CALL cpu_log( log_point(55), 'wind_turbine', 'stop' )
587
588       ENDIF       
589       
590!
[1]591!--    Start of intermediate step loop
592       intermediate_timestep_count = 0
593       DO  WHILE ( intermediate_timestep_count < &
594                   intermediate_timestep_count_max )
595
596          intermediate_timestep_count = intermediate_timestep_count + 1
597
598!
599!--       Set the steering factors for the prognostic equations which depend
600!--       on the timestep scheme
601          CALL timestep_scheme_steering
602
603!
[1128]604!--       Calculate those variables needed in the tendency terms which need
605!--       global communication
[1179]606          IF ( .NOT. use_single_reference_value  .AND. &
607               .NOT. use_initial_profile_as_reference )  THEN
608!
609!--          Horizontally averaged profiles to be used as reference state in
610!--          buoyancy terms (WARNING: only the respective last call of
611!--          calc_mean_profile defines the reference state!)
[1365]612             IF ( .NOT. neutral )  THEN
613                CALL calc_mean_profile( pt, 4 )
614                ref_state(:)  = hom(:,1,4,0) ! this is used in the buoyancy term
615             ENDIF
616             IF ( ocean )  THEN
[2031]617                CALL calc_mean_profile( rho_ocean, 64 )
[1365]618                ref_state(:)  = hom(:,1,64,0)
619             ENDIF
620             IF ( humidity )  THEN
621                CALL calc_mean_profile( vpt, 44 )
622                ref_state(:)  = hom(:,1,44,0)
623             ENDIF
[2617]624!
625!--          Assure that ref_state does not become zero at any level
626!--          ( might be the case if a vertical level is completely occupied
627!--            with topography ).
628             ref_state = MERGE( MAXVAL(ref_state), ref_state,                  &
629                                ref_state == 0.0_wp )
[1179]630          ENDIF
631
[1128]632          IF ( .NOT. constant_diffusion )  CALL production_e_init
633          IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
634               intermediate_timestep_count == 1 )  CALL ws_statistics
[1365]635!
636!--       In case of nudging calculate current nudging time scale and horizontal
[1380]637!--       means of u, v, pt and q
[1365]638          IF ( nudging )  THEN
639             CALL calc_tnudge( simulated_time )
640             CALL calc_mean_profile( u, 1 )
641             CALL calc_mean_profile( v, 2 )
642             CALL calc_mean_profile( pt, 4 )
643             CALL calc_mean_profile( q, 41 )
644          ENDIF
[1128]645
646!
[1]647!--       Solve the prognostic equations. A fast cache optimized version with
648!--       only one single loop is used in case of Piascek-Williams advection
649!--       scheme. NEC vector machines use a different version, because
650!--       in the other versions a good vectorization is prohibited due to
651!--       inlining problems.
[1019]652          IF ( loop_optimization == 'cache' )  THEN
653             CALL prognostic_equations_cache
654          ELSEIF ( loop_optimization == 'vector' )  THEN
[63]655             CALL prognostic_equations_vector
[1]656          ENDIF
657
658!
[849]659!--       Particle transport/physics with the Lagrangian particle model
660!--       (only once during intermediate steps, because it uses an Euler-step)
[1128]661!--       ### particle model should be moved before prognostic_equations, in order
662!--       to regard droplet interactions directly
[63]663          IF ( particle_advection  .AND.                         &
664               simulated_time >= particle_advection_start  .AND. &
[1]665               intermediate_timestep_count == 1 )  THEN
[849]666             CALL lpm
667             first_call_lpm = .FALSE.
[1]668          ENDIF
669
670!
671!--       Interaction of droplets with temperature and specific humidity.
672!--       Droplet condensation and evaporation is calculated within
673!--       advec_particles.
674          IF ( cloud_droplets  .AND.  &
675               intermediate_timestep_count == intermediate_timestep_count_max )&
676          THEN
677             CALL interaction_droplets_ptq
678          ENDIF
679
680!
681!--       Exchange of ghost points (lateral boundary conditions)
[2118]682          CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'start' )
[1113]683
[2118]684          CALL exchange_horiz( u_p, nbgp )
685          CALL exchange_horiz( v_p, nbgp )
686          CALL exchange_horiz( w_p, nbgp )
687          CALL exchange_horiz( pt_p, nbgp )
688          IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e_p, nbgp )
[2696]689          IF ( rans_tke_e  .OR.  wang_kernel  .OR.  collision_turbulence       &
690               .OR.  use_sgs_for_particles )  THEN
691             IF ( rans_tke_e )  THEN
692                CALL exchange_horiz( diss_p, nbgp )
693             ELSE
694                CALL exchange_horiz( diss, nbgp )
695             ENDIF
696          ENDIF
[2118]697          IF ( ocean )  THEN
698             CALL exchange_horiz( sa_p, nbgp )
699             CALL exchange_horiz( rho_ocean, nbgp )
700             CALL exchange_horiz( prho, nbgp )
701          ENDIF
702          IF ( humidity )  THEN
703             CALL exchange_horiz( q_p, nbgp )
[2292]704             IF ( cloud_physics .AND. microphysics_morrison )  THEN
705                CALL exchange_horiz( qc_p, nbgp )
706                CALL exchange_horiz( nc_p, nbgp )
707             ENDIF
[2118]708             IF ( cloud_physics .AND. microphysics_seifert )  THEN
709                CALL exchange_horiz( qr_p, nbgp )
710                CALL exchange_horiz( nr_p, nbgp )
[1053]711             ENDIF
[2118]712          ENDIF
713          IF ( cloud_droplets )  THEN
714             CALL exchange_horiz( ql, nbgp )
715             CALL exchange_horiz( ql_c, nbgp )
716             CALL exchange_horiz( ql_v, nbgp )
717             CALL exchange_horiz( ql_vp, nbgp )
718          ENDIF
[2696]719          IF ( passive_scalar )  CALL exchange_horiz( s_p, nbgp )
720          IF ( air_chemistry )  THEN
721             DO  n = 1, nspec     
722                CALL exchange_horiz( chem_species(n)%conc_p, nbgp ) 
723             ENDDO
[2118]724          ENDIF
[2766]725
[2118]726          CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
[1128]727
[1]728!
729!--       Boundary conditions for the prognostic quantities (except of the
730!--       velocities at the outflow in case of a non-cyclic lateral wall)
[1113]731          CALL boundary_conds
[1]732!
[73]733!--       Swap the time levels in preparation for the next time step.
734          CALL swap_timelevel
735
[2365]736!
737!--       Vertical nesting: Interpolate fine grid data to the coarse grid
738          IF ( vnest_init ) THEN
739             CALL cpu_log( log_point(81), 'vnest_anterpolate', 'start' )
740             CALL vnest_anterpolate
741             CALL cpu_log( log_point(81), 'vnest_anterpolate', 'stop' )
742          ENDIF
743
[1764]744          IF ( nested_run )  THEN
[1797]745
[1764]746             CALL cpu_log( log_point(60), 'nesting', 'start' )
[1762]747!
[1933]748!--          Domain nesting. The data transfer subroutines pmci_parent_datatrans
749!--          and pmci_child_datatrans are called inside the wrapper
[1797]750!--          subroutine pmci_datatrans according to the control parameters
751!--          nesting_mode and nesting_datatransfer_mode.
752!--          TO_DO: why is nesting_mode given as a parameter here?
753             CALL pmci_datatrans( nesting_mode )
[1762]754
[2174]755             IF ( TRIM( nesting_mode ) == 'two-way' .OR.                       &
[1933]756                  nesting_mode == 'vertical' )  THEN
[1762]757!
[1933]758!--             Exchange_horiz is needed for all parent-domains after the
[1764]759!--             anterpolation
760                CALL exchange_horiz( u, nbgp )
761                CALL exchange_horiz( v, nbgp )
762                CALL exchange_horiz( w, nbgp )
[2174]763                IF ( .NOT. neutral )  CALL exchange_horiz( pt, nbgp )
764
765                IF ( humidity )  THEN
766
767                   CALL exchange_horiz( q, nbgp )
768
[2292]769                   IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
770                       CALL exchange_horiz( qc, nbgp )
771                       CALL exchange_horiz( nc, nbgp )
772                   ENDIF
[2174]773                   IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
774                       CALL exchange_horiz( qr, nbgp )
775                       CALL exchange_horiz( nr, nbgp )
776                   ENDIF
777
778                ENDIF
779
780                IF ( passive_scalar )  CALL exchange_horiz( s, nbgp )
781                IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e, nbgp )
[2773]782
783                IF ( air_chemistry )  THEN
784                   DO  n = 1, nspec     
785                      CALL exchange_horiz( chem_species(n)%conc, nbgp ) 
786                   ENDDO
787                ENDIF
788
[1762]789             ENDIF
790!
[2311]791!--          Set boundary conditions again after interpolation and anterpolation.
792             CALL pmci_boundary_conds
793!
[1764]794!--          Correct the w top-BC in nest domains to ensure mass conservation.
[1933]795!--          This action must never be done for the root domain. Vertical
796!--          nesting implies mass conservation.
[1764]797             IF ( nest_domain )  THEN
798                CALL pmci_ensure_nest_mass_conservation
799             ENDIF
800
801             CALL cpu_log( log_point(60), 'nesting', 'stop' )
802
[1762]803          ENDIF
804
805!
[1]806!--       Temperature offset must be imposed at cyclic boundaries in x-direction
807!--       when a sloping surface is used
808          IF ( sloping_surface )  THEN
[707]809             IF ( nxl ==  0 )  pt(:,:,nxlg:nxl-1) = pt(:,:,nxlg:nxl-1) - &
810                                                    pt_slope_offset
811             IF ( nxr == nx )  pt(:,:,nxr+1:nxrg) = pt(:,:,nxr+1:nxrg) + &
812                                                    pt_slope_offset
[1]813          ENDIF
814
815!
[151]816!--       Impose a turbulent inflow using the recycling method
817          IF ( turbulent_inflow )  CALL  inflow_turbulence
818
819!
[2259]820!--       Impose a turbulent inflow using synthetic generated turbulence
[2776]821         IF ( use_syn_turb_gen ) THEN
[2259]822            CALL  stg_main
823         ENDIF
824
825!
[2050]826!--       Set values at outflow boundary using the special outflow condition
827          IF ( turbulent_outflow )  CALL  outflow_turbulence
828
829!
[1]830!--       Impose a random perturbation on the horizontal velocity field
[106]831          IF ( create_disturbances  .AND.                                      &
832               ( call_psolver_at_all_substeps  .AND.                           &
[1]833               intermediate_timestep_count == intermediate_timestep_count_max )&
[106]834          .OR. ( .NOT. call_psolver_at_all_substeps  .AND.                     &
835               intermediate_timestep_count == 1 ) )                            &
[1]836          THEN
837             time_disturb = time_disturb + dt_3d
838             IF ( time_disturb >= dt_disturb )  THEN
[1736]839                IF ( disturbance_energy_limit /= 0.0_wp  .AND.                 &
840                     hom(nzb+5,1,pr_palm,0) < disturbance_energy_limit )  THEN
[2232]841                   CALL disturb_field( 'u', tend, u )
842                   CALL disturb_field( 'v', tend, v )
[2178]843                ELSEIF ( ( .NOT. bc_lr_cyc  .OR.  .NOT. bc_ns_cyc )            &
[2696]844                     .AND. .NOT. nest_domain  .AND.  .NOT.  forcing )  THEN
[1]845!
846!--                Runs with a non-cyclic lateral wall need perturbations
847!--                near the inflow throughout the whole simulation
848                   dist_range = 1
[2232]849                   CALL disturb_field( 'u', tend, u )
850                   CALL disturb_field( 'v', tend, v )
[1]851                   dist_range = 0
852                ENDIF
853                time_disturb = time_disturb - dt_disturb
854             ENDIF
855          ENDIF
856
857!
[2696]858!--       Map forcing data derived from larger scale model onto domain
859!--       boundaries.
860          IF ( forcing  .AND.  intermediate_timestep_count ==                  &
861                               intermediate_timestep_count_max  )  THEN
862             CALL forcing_bc
863          ENDIF
864!
865!--       Moreover, ensure mass conservation at each RK substep.
866          IF ( forcing )  CALL forcing_bc_mass_conservation
867
868!
[1]869!--       Reduce the velocity divergence via the equation for perturbation
870!--       pressure.
[106]871          IF ( intermediate_timestep_count == 1  .OR. &
872                call_psolver_at_all_substeps )  THEN
[2365]873
874             IF (  vnest_init ) THEN
875!
876!--             Compute pressure in the CG, interpolate top boundary conditions
877!--             to the FG and then compute pressure in the FG
878                IF ( coupling_mode == 'vnested_crse' )  CALL pres
879
880                CALL cpu_log( log_point(82), 'vnest_bc', 'start' )
881                CALL vnest_boundary_conds
882                CALL cpu_log( log_point(82), 'vnest_bc', 'stop' )
883 
884                IF ( coupling_mode == 'vnested_fine' )  CALL pres
885
886!--             Anterpolate TKE, satisfy Germano Identity
887                CALL cpu_log( log_point(83), 'vnest_anter_e', 'start' )
888                CALL vnest_anterpolate_e
889                CALL cpu_log( log_point(83), 'vnest_anter_e', 'stop' )
890
891             ELSE
892
893                CALL pres
894
895             ENDIF
896
[1]897          ENDIF
898
899!
900!--       If required, compute liquid water content
[1015]901          IF ( cloud_physics )  THEN
902             CALL calc_liquid_water_content
903          ENDIF
[2174]904!
[1115]905!--       If required, compute virtual potential temperature
906          IF ( humidity )  THEN
907             CALL compute_vpt 
908          ENDIF 
[1585]909
[1]910!
911!--       Compute the diffusion quantities
912          IF ( .NOT. constant_diffusion )  THEN
913
914!
[1276]915!--          Determine surface fluxes shf and qsws and surface values
916!--          pt_surface and q_surface in dependence on data from external
917!--          file LSF_DATA respectively
918             IF ( ( large_scale_forcing .AND. lsf_surf ) .AND. &
919                 intermediate_timestep_count == intermediate_timestep_count_max )&
920             THEN
[2320]921                CALL ls_forcing_surf( simulated_time )
[1276]922             ENDIF
923
924!
[2232]925!--          First the vertical (and horizontal) fluxes in the surface
926!--          (constant flux) layer are computed
[1691]927             IF ( constant_flux_layer )  THEN
928                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'start' )
929                CALL surface_layer_fluxes
930                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'stop' )
[1]931             ENDIF
932!
[1691]933!--          If required, solve the energy balance for the surface and run soil
[2232]934!--          model. Call for horizontal as well as vertical surfaces
[2696]935             IF ( land_surface .AND. time_since_reference_point >= skip_time_do_lsm)  THEN
[1691]936
937                CALL cpu_log( log_point(54), 'land_surface', 'start' )
[2232]938!
939!--             Call for horizontal upward-facing surfaces
940                CALL lsm_energy_balance( .TRUE., -1 )
[2299]941                CALL lsm_soil_model( .TRUE., -1, .TRUE. )
[2232]942!
943!--             Call for northward-facing surfaces
944                CALL lsm_energy_balance( .FALSE., 0 )
[2299]945                CALL lsm_soil_model( .FALSE., 0, .TRUE. )
[2232]946!
947!--             Call for southward-facing surfaces
948                CALL lsm_energy_balance( .FALSE., 1 )
[2299]949                CALL lsm_soil_model( .FALSE., 1, .TRUE. )
[2232]950!
951!--             Call for eastward-facing surfaces
952                CALL lsm_energy_balance( .FALSE., 2 )
[2299]953                CALL lsm_soil_model( .FALSE., 2, .TRUE. )
[2232]954!
955!--             Call for westward-facing surfaces
956                CALL lsm_energy_balance( .FALSE., 3 )
[2299]957                CALL lsm_soil_model( .FALSE., 3, .TRUE. )
[2696]958!
959!--             At the end, set boundary conditons for potential temperature
960!--             and humidity after running the land-surface model. This
961!--             might be important for the nesting, where arrays are transfered.
962                CALL lsm_boundary_condition
[2232]963
[1691]964                CALL cpu_log( log_point(54), 'land_surface', 'stop' )
965             ENDIF
966!
[2007]967!--          If required, solve the energy balance for urban surfaces and run
968!--          the material heat model
969             IF (urban_surface) THEN
970                CALL cpu_log( log_point(74), 'urban_surface', 'start' )
971                CALL usm_surface_energy_balance
972                IF ( usm_material_model )  THEN
[2696]973                   CALL usm_green_heat_model
[2007]974                   CALL usm_material_heat_model
975                ENDIF
[2696]976
977                CALL usm_temperature_near_surface
978!
979!--             At the end, set boundary conditons for potential temperature
980!--             and humidity after running the urban-surface model. This
981!--             might be important for the nesting, where arrays are transfered.
982                CALL usm_boundary_condition
983
[2007]984                CALL cpu_log( log_point(74), 'urban_surface', 'stop' )
985             ENDIF
986!
[1]987!--          Compute the diffusion coefficients
988             CALL cpu_log( log_point(17), 'diffusivities', 'start' )
[75]989             IF ( .NOT. humidity ) THEN
[97]990                IF ( ocean )  THEN
[2696]991                   CALL tcm_diffusivities( prho, prho_reference )
[97]992                ELSE
[2696]993                   CALL tcm_diffusivities( pt, pt_reference )
[97]994                ENDIF
[1]995             ELSE
[2696]996                CALL tcm_diffusivities( vpt, pt_reference )
[1]997             ENDIF
998             CALL cpu_log( log_point(17), 'diffusivities', 'stop' )
[2696]999!
1000!--          Vertical nesting: set fine grid eddy viscosity top boundary condition
1001             IF ( vnest_init )  CALL vnest_boundary_conds_khkm
[1]1002
1003          ENDIF
1004
[1691]1005!
1006!--       If required, calculate radiative fluxes and heating rates
1007          IF ( radiation .AND. intermediate_timestep_count                     &
[2299]1008               == intermediate_timestep_count_max .AND. time_since_reference_point >    &
[1691]1009               skip_time_do_radiation )  THEN
1010
1011               time_radiation = time_radiation + dt_3d
1012
1013             IF ( time_radiation >= dt_radiation .OR. force_radiation_call )   &
1014             THEN
1015
1016                CALL cpu_log( log_point(50), 'radiation', 'start' )
1017
1018                IF ( .NOT. force_radiation_call )  THEN
1019                   time_radiation = time_radiation - dt_radiation
1020                ENDIF
1021
[1976]1022                CALL radiation_control
[1691]1023
1024                CALL cpu_log( log_point(50), 'radiation', 'stop' )
[2007]1025
[2696]1026                IF ( urban_surface  .OR.  land_surface  .AND.                  &
1027                     radiation_interactions )  THEN
1028                   CALL cpu_log( log_point(75), 'radiation_interaction', 'start' )
1029                   CALL radiation_interaction
1030                   CALL cpu_log( log_point(75), 'radiation_interaction', 'stop' )
[2007]1031                ENDIF
1032
[1691]1033             ENDIF
1034          ENDIF
1035
[1]1036       ENDDO   ! Intermediate step loop
1037!
[2766]1038!--    If required, consider chemical emissions
1039!--    (todo (FK): Implement hourly call of emissions, using time_utc from
1040!--                data_and_time_mod.f90;
1041!--                move the CALL to appropriate location)
1042       IF ( air_chemistry ) THEN
1043          CALL chem_emissions
1044       ENDIF
[2696]1045!
1046!--    If required, do UV exposure calculations
1047       IF ( uv_exposure )  THEN
1048          CALL uvem_calc_exposure
1049       ENDIF
1050!
[1]1051!--    Increase simulation time and output times
[1111]1052       nr_timesteps_this_run      = nr_timesteps_this_run + 1
[291]1053       current_timestep_number    = current_timestep_number + 1
1054       simulated_time             = simulated_time   + dt_3d
1055       simulated_time_chr         = time_to_string( simulated_time )
1056       time_since_reference_point = simulated_time - coupling_start_time
1057
[1957]1058
1059
[1]1060       IF ( simulated_time >= skip_time_data_output_av )  THEN
1061          time_do_av         = time_do_av       + dt_3d
1062       ENDIF
1063       IF ( simulated_time >= skip_time_do2d_xy )  THEN
1064          time_do2d_xy       = time_do2d_xy     + dt_3d
1065       ENDIF
1066       IF ( simulated_time >= skip_time_do2d_xz )  THEN
1067          time_do2d_xz       = time_do2d_xz     + dt_3d
1068       ENDIF
1069       IF ( simulated_time >= skip_time_do2d_yz )  THEN
1070          time_do2d_yz       = time_do2d_yz     + dt_3d
1071       ENDIF
1072       IF ( simulated_time >= skip_time_do3d    )  THEN
1073          time_do3d          = time_do3d        + dt_3d
1074       ENDIF
[410]1075       DO  mid = 1, masks
1076          IF ( simulated_time >= skip_time_domask(mid) )  THEN
1077             time_domask(mid)= time_domask(mid) + dt_3d
1078          ENDIF
1079       ENDDO
[1]1080       time_dvrp          = time_dvrp        + dt_3d
1081       IF ( simulated_time >= skip_time_dosp )  THEN
1082          time_dosp       = time_dosp        + dt_3d
1083       ENDIF
1084       time_dots          = time_dots        + dt_3d
[849]1085       IF ( .NOT. first_call_lpm )  THEN
[1]1086          time_dopts      = time_dopts       + dt_3d
1087       ENDIF
1088       IF ( simulated_time >= skip_time_dopr )  THEN
1089          time_dopr       = time_dopr        + dt_3d
1090       ENDIF
1091       time_dopr_listing          = time_dopr_listing        + dt_3d
1092       time_run_control   = time_run_control + dt_3d
1093
1094!
[102]1095!--    Data exchange between coupled models
[2365]1096       IF ( coupling_mode /= 'uncoupled'  .AND.  run_coupled                   &
1097                                          .AND. .NOT. vnested )  THEN
[102]1098          time_coupling = time_coupling + dt_3d
[343]1099
[108]1100!
1101!--       In case of model termination initiated by the local model
1102!--       (terminate_coupled > 0), the coupler must be skipped because it would
1103!--       cause an MPI intercomminucation hang.
1104!--       If necessary, the coupler will be called at the beginning of the
1105!--       next restart run.
1106          DO WHILE ( time_coupling >= dt_coupling .AND. terminate_coupled == 0 )
[102]1107             CALL surface_coupler
1108             time_coupling = time_coupling - dt_coupling
1109          ENDDO
1110       ENDIF
1111
1112!
[2817]1113!--    Execute the gust module actions
1114       IF ( gust_module_enabled )  THEN
1115          CALL gust_actions( 'after_integration' )
1116       ENDIF
1117
1118!
[46]1119!--    Execute user-defined actions
1120       CALL user_actions( 'after_integration' )
1121
1122!
[1]1123!--    If Galilei transformation is used, determine the distance that the
1124!--    model has moved so far
1125       IF ( galilei_transformation )  THEN
1126          advected_distance_x = advected_distance_x + u_gtrans * dt_3d
1127          advected_distance_y = advected_distance_y + v_gtrans * dt_3d
1128       ENDIF
1129
1130!
1131!--    Check, if restart is necessary (because cpu-time is expiring or
1132!--    because it is forced by user) and set stop flag
[108]1133!--    This call is skipped if the remote model has already initiated a restart.
1134       IF ( .NOT. terminate_run )  CALL check_for_restart
[1]1135
1136!
1137!--    Carry out statistical analysis and output at the requested output times.
1138!--    The MOD function is used for calculating the output time counters (like
1139!--    time_dopr) in order to regard a possible decrease of the output time
1140!--    interval in case of restart runs
1141
1142!
1143!--    Set a flag indicating that so far no statistics have been created
1144!--    for this time step
1145       flow_statistics_called = .FALSE.
1146
1147!
1148!--    If required, call flow_statistics for averaging in time
[1342]1149       IF ( averaging_interval_pr /= 0.0_wp  .AND.  &
[1]1150            ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.  &
1151            simulated_time >= skip_time_dopr )  THEN
1152          time_dopr_av = time_dopr_av + dt_3d
1153          IF ( time_dopr_av >= dt_averaging_input_pr )  THEN
1154             do_sum = .TRUE.
1155             time_dopr_av = MOD( time_dopr_av, &
1156                                    MAX( dt_averaging_input_pr, dt_3d ) )
1157          ENDIF
1158       ENDIF
1159       IF ( do_sum )  CALL flow_statistics
1160
1161!
[410]1162!--    Sum-up 3d-arrays for later output of time-averaged 2d/3d/masked data
[1342]1163       IF ( averaging_interval /= 0.0_wp  .AND.                                &
[1]1164            ( dt_data_output_av - time_do_av ) <= averaging_interval  .AND. &
1165            simulated_time >= skip_time_data_output_av )                    &
1166       THEN
1167          time_do_sla = time_do_sla + dt_3d
1168          IF ( time_do_sla >= dt_averaging_input )  THEN
1169             CALL sum_up_3d_data
1170             average_count_3d = average_count_3d + 1
1171             time_do_sla = MOD( time_do_sla, MAX( dt_averaging_input, dt_3d ) )
1172          ENDIF
1173       ENDIF
1174
1175!
1176!--    Calculate spectra for time averaging
[1342]1177       IF ( averaging_interval_sp /= 0.0_wp  .AND.  &
[1]1178            ( dt_dosp - time_dosp ) <= averaging_interval_sp  .AND.  &
1179            simulated_time >= skip_time_dosp )  THEN
1180          time_dosp_av = time_dosp_av + dt_3d
1181          IF ( time_dosp_av >= dt_averaging_input_pr )  THEN
1182             CALL calc_spectra
1183             time_dosp_av = MOD( time_dosp_av, &
1184                                 MAX( dt_averaging_input_pr, dt_3d ) )
1185          ENDIF
1186       ENDIF
1187
1188!
[1957]1189!--    Call flight module and output data
1190       IF ( virtual_flight )  THEN
1191          CALL flight_measurement
1192          CALL data_output_flight
1193       ENDIF
1194
1195!
[1]1196!--    Profile output (ASCII) on file
1197       IF ( time_dopr_listing >= dt_dopr_listing )  THEN
1198          CALL print_1d
1199          time_dopr_listing = MOD( time_dopr_listing, MAX( dt_dopr_listing, &
1200                                                           dt_3d ) )
1201       ENDIF
1202
1203!
1204!--    Graphic output for PROFIL
1205       IF ( time_dopr >= dt_dopr )  THEN
1206          IF ( dopr_n /= 0 )  CALL data_output_profiles
1207          time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
[1342]1208          time_dopr_av = 0.0_wp    ! due to averaging (see above)
[1]1209       ENDIF
1210
1211!
1212!--    Graphic output for time series
1213       IF ( time_dots >= dt_dots )  THEN
[48]1214          CALL data_output_tseries
[1]1215          time_dots = MOD( time_dots, MAX( dt_dots, dt_3d ) )
1216       ENDIF
1217
1218!
1219!--    Output of spectra (formatted for use with PROFIL), in case of no
1220!--    time averaging, spectra has to be calculated before
1221       IF ( time_dosp >= dt_dosp )  THEN
1222          IF ( average_count_sp == 0 )  CALL calc_spectra
1223          CALL data_output_spectra
1224          time_dosp = MOD( time_dosp, MAX( dt_dosp, dt_3d ) )
1225       ENDIF
1226
1227!
1228!--    2d-data output (cross-sections)
1229       IF ( time_do2d_xy >= dt_do2d_xy )  THEN
1230          CALL data_output_2d( 'xy', 0 )
1231          time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
1232       ENDIF
1233       IF ( time_do2d_xz >= dt_do2d_xz )  THEN
1234          CALL data_output_2d( 'xz', 0 )
1235          time_do2d_xz = MOD( time_do2d_xz, MAX( dt_do2d_xz, dt_3d ) )
1236       ENDIF
1237       IF ( time_do2d_yz >= dt_do2d_yz )  THEN
1238          CALL data_output_2d( 'yz', 0 )
1239          time_do2d_yz = MOD( time_do2d_yz, MAX( dt_do2d_yz, dt_3d ) )
1240       ENDIF
1241
1242!
1243!--    3d-data output (volume data)
1244       IF ( time_do3d >= dt_do3d )  THEN
1245          CALL data_output_3d( 0 )
1246          time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
1247       ENDIF
1248
1249!
[1783]1250!--    Masked data output
[410]1251       DO  mid = 1, masks
1252          IF ( time_domask(mid) >= dt_domask(mid) )  THEN
1253             CALL data_output_mask( 0 )
1254             time_domask(mid) = MOD( time_domask(mid),  &
1255                                     MAX( dt_domask(mid), dt_3d ) )
1256          ENDIF
1257       ENDDO
1258
1259!
1260!--    Output of time-averaged 2d/3d/masked data
[1]1261       IF ( time_do_av >= dt_data_output_av )  THEN
1262          CALL average_3d_data
1263          CALL data_output_2d( 'xy', 1 )
1264          CALL data_output_2d( 'xz', 1 )
1265          CALL data_output_2d( 'yz', 1 )
1266          CALL data_output_3d( 1 )
[410]1267          DO  mid = 1, masks
1268             CALL data_output_mask( 1 )
1269          ENDDO
[1]1270          time_do_av = MOD( time_do_av, MAX( dt_data_output_av, dt_3d ) )
1271       ENDIF
1272
1273!
1274!--    Output of particle time series
[253]1275       IF ( particle_advection )  THEN
1276          IF ( time_dopts >= dt_dopts  .OR. &
1277               ( simulated_time >= particle_advection_start  .AND. &
[849]1278                 first_call_lpm ) )  THEN
[253]1279             CALL data_output_ptseries
1280             time_dopts = MOD( time_dopts, MAX( dt_dopts, dt_3d ) )
1281          ENDIF
[1]1282       ENDIF
1283
1284!
1285!--    Output of dvrp-graphics (isosurface, particles, slicer)
1286#if defined( __dvrp_graphics )
1287       CALL DVRP_LOG_EVENT( -2, current_timestep_number-1 )
1288#endif
1289       IF ( time_dvrp >= dt_dvrp )  THEN
1290          CALL data_output_dvrp
1291          time_dvrp = MOD( time_dvrp, MAX( dt_dvrp, dt_3d ) )
1292       ENDIF
1293#if defined( __dvrp_graphics )
1294       CALL DVRP_LOG_EVENT( 2, current_timestep_number )
1295#endif
1296
1297!
1298!--    If required, set the heat flux for the next time step at a random value
[2232]1299       IF ( constant_heatflux  .AND.  random_heatflux )  THEN
1300          IF ( surf_def_h(0)%ns >= 1 )  CALL disturb_heatflux( surf_def_h(0) )
1301          IF ( surf_lsm_h%ns    >= 1 )  CALL disturb_heatflux( surf_lsm_h    )
1302          IF ( surf_usm_h%ns    >= 1 )  CALL disturb_heatflux( surf_usm_h    )
1303       ENDIF
[1]1304
1305!
[2817]1306!--    Execute the gust module actions
1307       IF ( gust_module_enabled )  THEN
1308          CALL gust_actions( 'after_timestep' )
1309       ENDIF
1310
1311!
[1]1312!--    Execute user-defined actions
1313       CALL user_actions( 'after_timestep' )
1314
[1402]1315!
[1918]1316!--    Determine size of next time step. Save timestep dt_3d because it is
1317!--    newly calculated in routine timestep, but required further below for
1318!--    steering the run control output interval
1319       dt_3d_old = dt_3d
1320       CALL timestep
1321
1322!
[1925]1323!--    Synchronize the timestep in case of nested run.
1324       IF ( nested_run )  THEN
1325!
1326!--       Synchronize by unifying the time step.
1327!--       Global minimum of all time-steps is used for all.
1328          CALL pmci_synchronize
1329       ENDIF
1330
1331!
[1918]1332!--    Computation and output of run control parameters.
1333!--    This is also done whenever perturbations have been imposed
1334       IF ( time_run_control >= dt_run_control  .OR.                     &
1335            timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created ) &
1336       THEN
1337          CALL run_control
1338          IF ( time_run_control >= dt_run_control )  THEN
1339             time_run_control = MOD( time_run_control, &
1340                                     MAX( dt_run_control, dt_3d_old ) )
1341          ENDIF
1342       ENDIF
1343
1344!
[1402]1345!--    Output elapsed simulated time in form of a progress bar on stdout
1346       IF ( myid == 0 )  CALL output_progress_bar
1347
[1]1348       CALL cpu_log( log_point_s(10), 'timesteps', 'stop' )
1349
[667]1350
[1]1351    ENDDO   ! time loop
1352
[2365]1353!-- Vertical nesting: Deallocate variables initialized for vertical nesting   
1354    IF ( vnest_init )  CALL vnest_deallocate
1355
[1402]1356    IF ( myid == 0 )  CALL finish_progress_bar
1357
[1]1358#if defined( __dvrp_graphics )
1359    CALL DVRP_LOG_EVENT( -2, current_timestep_number )
1360#endif
1361
[1402]1362    CALL location_message( 'finished time-stepping', .TRUE. )
[1384]1363
[1]1364 END SUBROUTINE time_integration
Note: See TracBrowser for help on using the repository browser.