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

Last change on this file since 2008 was 2008, checked in by kanani, 8 years ago

last commit documented

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