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

Last change on this file since 3494 was 3484, checked in by hellstea, 5 years ago

Nesting interpolation made mass-conservative

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