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

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

changes in the course of urban surface model implementation

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