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

Last change on this file since 1165 was 1132, checked in by raasch, 12 years ago

r1028 documented

  • Property svn:keywords set to Id
File size: 30.5 KB
Line 
1 SUBROUTINE time_integration
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later 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-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: time_integration.f90 1132 2013-04-12 14:35:30Z raasch $
27!
28! 1128 2013-04-12 06:19:32Z raasch
29! asynchronous transfer of ghost point data realized for acc-optimized version:
30! prognostic_equations are first called two times for those points required for
31! the left-right and north-south exchange, respectively, and then for the
32! remaining points,
33! those parts requiring global communication moved from prognostic_equations to
34! here
35!
36! 1115 2013-03-26 18:16:16Z hoffmann
37! calculation of qr and nr is restricted to precipitation
38!
39! 1113 2013-03-10 02:48:14Z raasch
40! GPU-porting of boundary conditions,
41! openACC directives updated
42! formal parameter removed from routine boundary_conds
43!
44! 1111 2013-03-08 23:54:10Z raasch
45! +internal timestep counter for cpu statistics added,
46! openACC directives updated
47!
48! 1092 2013-02-02 11:24:22Z raasch
49! unused variables removed
50!
51! 1065 2012-11-22 17:42:36Z hoffmann
52! exchange of diss (dissipation rate) in case of turbulence = .TRUE. added
53!
54! 1053 2012-11-13 17:11:03Z hoffmann
55! exchange of ghost points for nr, qr added
56!
57! 1036 2012-10-22 13:43:42Z raasch
58! code put under GPL (PALM 3.9)
59!
60! 1019 2012-09-28 06:46:45Z raasch
61! non-optimized version of prognostic_equations removed
62!
63! 1015 2012-09-27 09:23:24Z raasch
64! +call of prognostic_equations_acc
65!
66! 1001 2012-09-13 14:08:46Z raasch
67! all actions concerning leapfrog- and upstream-spline-scheme removed
68!
69! 849 2012-03-15 10:35:09Z raasch
70! advec_particles renamed lpm, first_call_advec_particles renamed first_call_lpm
71!
72! 825 2012-02-19 03:03:44Z raasch
73! wang_collision_kernel renamed wang_kernel
74!
75! 790 2011-11-29 03:11:20Z raasch
76! exchange of ghostpoints for array diss
77!
78! 707 2011-03-29 11:39:40Z raasch
79! bc_lr/ns replaced by bc_lr/ns_cyc, calls of exchange_horiz are modified,
80! adaption to sloping surface
81!
82! 667  2010-12-23 12:06:00Z suehring/gryschka
83! Calls of exchange_horiz are modified.
84! Adaption to slooping surface.
85!
86! 449 2010-02-02 11:23:59Z raasch
87! Bugfix: exchange of ghost points for prho included
88!
89! 410 2009-12-04 17:05:40Z letzel
90! masked data output
91!
92! 388 2009-09-23 09:40:33Z raasch
93! Using prho instead of rho in diffusvities.
94! Coupling with independent precursor runs.
95! Bugfix: output of particle time series only if particle advection is switched
96!         on
97!
98! 151 2008-03-07 13:42:18Z raasch
99! inflow turbulence is imposed by calling new routine inflow_turbulence
100!
101! 108 2007-08-24 15:10:38Z letzel
102! Call of new routine surface_coupler,
103! presure solver is called after the first Runge-Kutta substep instead of the
104! last in case that call_psolver_at_all_substeps = .F.; for this case, the
105! random perturbation has to be added to the velocity fields also after the
106! first substep
107!
108! 97 2007-06-21 08:23:15Z raasch
109! diffusivities is called with argument rho in case of ocean runs,
110! new argument pt_/prho_reference in calls of diffusivities,
111! ghostpoint exchange for salinity and density
112!
113! 87 2007-05-22 15:46:47Z raasch
114! var_hom renamed pr_palm
115!
116! 75 2007-03-22 09:54:05Z raasch
117! Move call of user_actions( 'after_integration' ) below increment of times
118! and counters,
119! calls of prognostic_equations_.. changed to .._noopt, .._cache, and
120! .._vector, these calls are now controlled by switch loop_optimization,
121! uxrp, vynp eliminated, 2nd+3rd argument removed from exchange horiz,
122! moisture renamed humidity
123!
124! RCS Log replace by Id keyword, revision history cleaned up
125!
126! Revision 1.8  2006/08/22 14:16:05  raasch
127! Disturbances are imposed only for the last Runge-Kutta-substep
128!
129! Revision 1.2  2004/04/30 13:03:40  raasch
130! decalpha-specific warning removed, routine name changed to time_integration,
131! particle advection is carried out only once during the intermediate steps,
132! impulse_advec renamed momentum_advec
133!
134! Revision 1.1  1997/08/11 06:19:04  raasch
135! Initial revision
136!
137!
138! Description:
139! ------------
140! Integration in time of the model equations, statistical analysis and graphic
141! output
142!------------------------------------------------------------------------------!
143
144    USE advec_ws
145    USE arrays_3d
146    USE averaging
147    USE buoyancy_mod
148    USE control_parameters
149    USE cpulog
150#if defined( __dvrp_graphics )
151    USE DVRP
152#endif
153    USE grid_variables
154    USE indices
155    USE interaction_droplets_ptq_mod
156    USE interfaces
157    USE particle_attributes
158    USE pegrid
159    USE production_e_mod
160    USE prognostic_equations_mod
161    USE statistics
162    USE user_actions_mod
163
164    IMPLICIT NONE
165
166    CHARACTER (LEN=9) ::  time_to_string
167
168!
169!-- At the beginning of a simulation determine the time step as well as
170!-- determine and print out the run control parameters
171    IF ( simulated_time == 0.0 )  CALL timestep
172
173    CALL run_control
174
175
176!
177!-- Data exchange between coupled models in case that a call has been omitted
178!-- at the end of the previous run of a job chain.
179    IF ( coupling_mode /= 'uncoupled'  .AND.  run_coupled )  THEN
180!
181!--    In case of model termination initiated by the local model the coupler
182!--    must not be called because this would again cause an MPI hang.
183       DO WHILE ( time_coupling >= dt_coupling .AND. terminate_coupled == 0 )
184          CALL surface_coupler
185          time_coupling = time_coupling - dt_coupling
186       ENDDO
187       IF (time_coupling == 0.0 .AND. time_since_reference_point < dt_coupling)&
188       THEN
189          time_coupling = time_since_reference_point
190       ENDIF
191    ENDIF
192
193
194#if defined( __dvrp_graphics )
195!
196!-- Time measurement with dvrp software 
197    CALL DVRP_LOG_EVENT( 2, current_timestep_number )
198#endif
199
200!
201!-- Start of the time loop
202    DO  WHILE ( simulated_time < end_time  .AND.  .NOT. stop_dt  .AND. &
203                .NOT. terminate_run )
204
205       CALL cpu_log( log_point_s(10), 'timesteps', 'start' )
206!
207!--    Determine size of next time step
208       IF ( simulated_time /= 0.0 )  CALL timestep
209!
210!--    Execute the user-defined actions
211       CALL user_actions( 'before_timestep' )
212
213!
214!--    Start of intermediate step loop
215       intermediate_timestep_count = 0
216       DO  WHILE ( intermediate_timestep_count < &
217                   intermediate_timestep_count_max )
218
219          intermediate_timestep_count = intermediate_timestep_count + 1
220
221!
222!--       Set the steering factors for the prognostic equations which depend
223!--       on the timestep scheme
224          CALL timestep_scheme_steering
225
226!
227!--       Calculate those variables needed in the tendency terms which need
228!--       global communication
229          IF ( .NOT. neutral )  CALL calc_mean_profile( pt, 4 )
230          IF ( ocean         )  CALL calc_mean_profile( rho, 64 )
231          IF ( humidity      )  CALL calc_mean_profile( vpt, 44 )
232          IF ( .NOT. constant_diffusion )  CALL production_e_init
233          IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
234               intermediate_timestep_count == 1 )  CALL ws_statistics
235
236!
237!--       Solve the prognostic equations. A fast cache optimized version with
238!--       only one single loop is used in case of Piascek-Williams advection
239!--       scheme. NEC vector machines use a different version, because
240!--       in the other versions a good vectorization is prohibited due to
241!--       inlining problems.
242          IF ( loop_optimization == 'cache' )  THEN
243             CALL prognostic_equations_cache
244          ELSEIF ( loop_optimization == 'vector' )  THEN
245             CALL prognostic_equations_vector
246          ELSEIF ( loop_optimization == 'acc' )  THEN
247!             i_left  = nxl;         i_right = nxr
248!             j_south = nys;         j_north = nyn
249!             CALL prognostic_equations_acc
250
251             i_left  = nxl;         i_right = nxl+nbgp-1
252             j_south = nys;         j_north = nyn
253             CALL prognostic_equations_acc
254             i_left  = nxr-nbgp+1;  i_right = nxr
255             j_south = nys;         j_north = nyn
256             CALL prognostic_equations_acc
257
258!
259!--          Exchange of ghost points (lateral boundary conditions)
260             IF ( background_communication )  THEN
261
262                CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'start' )
263               
264                send_receive = 'lr'
265                sendrecv_in_background = .TRUE.
266                req          = 0
267                req_count    = 0
268
269                IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
270                   on_device = .TRUE.         ! to be removed after complete porting
271                ELSE                          ! of ghost point exchange
272                   !$acc update host( e_p, pt_p, u_p, v_p, w_p )
273                ENDIF
274
275                CALL exchange_horiz( u_p, nbgp )
276                CALL exchange_horiz( v_p, nbgp )
277                CALL exchange_horiz( w_p, nbgp )
278                CALL exchange_horiz( pt_p, nbgp )
279                IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e_p, nbgp )
280                IF ( ocean )  THEN
281                   CALL exchange_horiz( sa_p, nbgp )
282                   CALL exchange_horiz( rho, nbgp )
283                  CALL exchange_horiz( prho, nbgp )
284                ENDIF
285                IF (humidity  .OR.  passive_scalar)  THEN
286                   CALL exchange_horiz( q_p, nbgp )
287                   IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
288                      CALL exchange_horiz( qr_p, nbgp )
289                      CALL exchange_horiz( nr_p, nbgp )
290                   ENDIF
291                ENDIF
292                IF ( cloud_droplets )  THEN
293                   CALL exchange_horiz( ql, nbgp )
294                   CALL exchange_horiz( ql_c, nbgp )
295                   CALL exchange_horiz( ql_v, nbgp )
296                   CALL exchange_horiz( ql_vp, nbgp )
297                ENDIF
298                IF ( wang_kernel  .OR.  turbulence )  CALL exchange_horiz( diss, nbgp )
299
300                IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
301                   on_device = .FALSE.        ! to be removed after complete porting
302                ELSE                          ! of ghost point exchange
303                   !$acc update device( e_p, pt_p, u_p, v_p, w_p )
304                ENDIF
305
306                sendrecv_in_background = .FALSE.
307
308                CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'pause' )
309
310             ENDIF
311
312             i_left  = nxl+nbgp;    i_right = nxr-nbgp
313             j_south = nys;         j_north = nys+nbgp-1
314             CALL prognostic_equations_acc
315             i_left  = nxl+nbgp;    i_right = nxr-nbgp
316             j_south = nyn-nbgp+1;  j_north = nyn
317             CALL prognostic_equations_acc
318
319             IF ( background_communication )  THEN
320                CALL cpu_log( log_point(41), 'exchange-horiz-wait', 'start' )
321#if defined( __parallel )
322                CALL MPI_WAITALL( req_count, req, wait_stat, ierr )
323#endif
324                CALL cpu_log( log_point(41), 'exchange-horiz-wait', 'pause' )
325
326                CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'continue' )
327
328                send_receive = 'ns'
329                sendrecv_in_background = .TRUE.
330                req          = 0
331                req_count    = 0
332
333                IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
334                   on_device = .TRUE.         ! to be removed after complete porting
335                ELSE                          ! of ghost point exchange
336                   !$acc update host( e_p, pt_p, u_p, v_p, w_p )
337                ENDIF
338
339                CALL exchange_horiz( u_p, nbgp )
340                CALL exchange_horiz( v_p, nbgp )
341                CALL exchange_horiz( w_p, nbgp )
342                CALL exchange_horiz( pt_p, nbgp )
343                IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e_p, nbgp )
344                IF ( ocean )  THEN
345                   CALL exchange_horiz( sa_p, nbgp )
346                   CALL exchange_horiz( rho, nbgp )
347                  CALL exchange_horiz( prho, nbgp )
348                ENDIF
349                IF (humidity  .OR.  passive_scalar)  THEN
350                   CALL exchange_horiz( q_p, nbgp )
351                   IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
352                      CALL exchange_horiz( qr_p, nbgp )
353                      CALL exchange_horiz( nr_p, nbgp )
354                   ENDIF
355                ENDIF
356                IF ( cloud_droplets )  THEN
357                   CALL exchange_horiz( ql, nbgp )
358                   CALL exchange_horiz( ql_c, nbgp )
359                   CALL exchange_horiz( ql_v, nbgp )
360                   CALL exchange_horiz( ql_vp, nbgp )
361                ENDIF
362                IF ( wang_kernel  .OR.  turbulence )  CALL exchange_horiz( diss, nbgp )
363
364                IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
365                   on_device = .FALSE.        ! to be removed after complete porting
366                ELSE                          ! of ghost point exchange
367                   !$acc update device( e_p, pt_p, u_p, v_p, w_p )
368                ENDIF
369
370                sendrecv_in_background = .FALSE.
371
372                CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
373
374             ENDIF
375
376             i_left  = nxl+nbgp;    i_right = nxr-nbgp
377             j_south = nys+nbgp;    j_north = nyn-nbgp
378             CALL prognostic_equations_acc
379
380             IF ( background_communication )  THEN
381                CALL cpu_log( log_point(41), 'exchange-horiz-wait', 'continue' )
382#if defined( __parallel )
383                CALL MPI_WAITALL( req_count, req, wait_stat, ierr )
384#endif
385                send_receive = 'al'
386                CALL cpu_log( log_point(41), 'exchange-horiz-wait', 'stop' )
387             ENDIF
388
389          ENDIF
390
391!
392!--       Particle transport/physics with the Lagrangian particle model
393!--       (only once during intermediate steps, because it uses an Euler-step)
394!--       ### particle model should be moved before prognostic_equations, in order
395!--       to regard droplet interactions directly
396          IF ( particle_advection  .AND.                         &
397               simulated_time >= particle_advection_start  .AND. &
398               intermediate_timestep_count == 1 )  THEN
399             CALL lpm
400             first_call_lpm = .FALSE.
401          ENDIF
402
403!
404!--       Interaction of droplets with temperature and specific humidity.
405!--       Droplet condensation and evaporation is calculated within
406!--       advec_particles.
407          IF ( cloud_droplets  .AND.  &
408               intermediate_timestep_count == intermediate_timestep_count_max )&
409          THEN
410             CALL interaction_droplets_ptq
411          ENDIF
412
413!
414!--       Exchange of ghost points (lateral boundary conditions)
415          IF ( .NOT. background_communication )  THEN
416
417             CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'start' )
418
419             IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
420                on_device = .TRUE.         ! to be removed after complete porting
421             ELSE                          ! of ghost point exchange
422                !$acc update host( e_p, pt_p, u_p, v_p, w_p )
423             ENDIF
424
425             CALL exchange_horiz( u_p, nbgp )
426             CALL exchange_horiz( v_p, nbgp )
427             CALL exchange_horiz( w_p, nbgp )
428             CALL exchange_horiz( pt_p, nbgp )
429             IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e_p, nbgp )
430             IF ( ocean )  THEN
431                CALL exchange_horiz( sa_p, nbgp )
432                CALL exchange_horiz( rho, nbgp )
433                CALL exchange_horiz( prho, nbgp )
434             ENDIF
435             IF (humidity  .OR.  passive_scalar)  THEN
436                CALL exchange_horiz( q_p, nbgp )
437                IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
438                   CALL exchange_horiz( qr_p, nbgp )
439                   CALL exchange_horiz( nr_p, nbgp )
440                ENDIF
441             ENDIF
442             IF ( cloud_droplets )  THEN
443                CALL exchange_horiz( ql, nbgp )
444                CALL exchange_horiz( ql_c, nbgp )
445                CALL exchange_horiz( ql_v, nbgp )
446                CALL exchange_horiz( ql_vp, nbgp )
447             ENDIF
448             IF ( wang_kernel  .OR.  turbulence )  CALL exchange_horiz( diss, nbgp )
449
450             IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
451                on_device = .FALSE.        ! to be removed after complete porting
452             ELSE                          ! of ghost point exchange
453                !$acc update device( e_p, pt_p, u_p, v_p, w_p )
454             ENDIF
455
456             CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
457
458          ENDIF
459
460!
461!--       Boundary conditions for the prognostic quantities (except of the
462!--       velocities at the outflow in case of a non-cyclic lateral wall)
463          CALL boundary_conds
464
465!
466!--       Swap the time levels in preparation for the next time step.
467          CALL swap_timelevel
468
469!
470!--       Temperature offset must be imposed at cyclic boundaries in x-direction
471!--       when a sloping surface is used
472          IF ( sloping_surface )  THEN
473             IF ( nxl ==  0 )  pt(:,:,nxlg:nxl-1) = pt(:,:,nxlg:nxl-1) - &
474                                                    pt_slope_offset
475             IF ( nxr == nx )  pt(:,:,nxr+1:nxrg) = pt(:,:,nxr+1:nxrg) + &
476                                                    pt_slope_offset
477          ENDIF
478
479!
480!--       Impose a turbulent inflow using the recycling method
481          IF ( turbulent_inflow )  CALL  inflow_turbulence
482
483!
484!--       Impose a random perturbation on the horizontal velocity field
485          IF ( create_disturbances  .AND.                                      &
486               ( call_psolver_at_all_substeps  .AND.                           &
487               intermediate_timestep_count == intermediate_timestep_count_max )&
488          .OR. ( .NOT. call_psolver_at_all_substeps  .AND.                     &
489               intermediate_timestep_count == 1 ) )                            &
490          THEN
491             time_disturb = time_disturb + dt_3d
492             IF ( time_disturb >= dt_disturb )  THEN
493                !$acc update host( u, v )
494                IF ( numprocs == 1 )  on_device = .FALSE.  ! workaround, remove later
495                IF ( hom(nzb+5,1,pr_palm,0) < disturbance_energy_limit )  THEN
496                   CALL disturb_field( nzb_u_inner, tend, u )
497                   CALL disturb_field( nzb_v_inner, tend, v )
498                ELSEIF ( .NOT. bc_lr_cyc  .OR.  .NOT. bc_ns_cyc )  THEN
499!
500!--                Runs with a non-cyclic lateral wall need perturbations
501!--                near the inflow throughout the whole simulation
502                   dist_range = 1
503                   CALL disturb_field( nzb_u_inner, tend, u )
504                   CALL disturb_field( nzb_v_inner, tend, v )
505                   dist_range = 0
506                ENDIF
507                IF ( numprocs == 1 )  on_device = .TRUE.  ! workaround, remove later
508                !$acc update device( u, v )
509                time_disturb = time_disturb - dt_disturb
510             ENDIF
511          ENDIF
512
513!
514!--       Reduce the velocity divergence via the equation for perturbation
515!--       pressure.
516          IF ( intermediate_timestep_count == 1  .OR. &
517                call_psolver_at_all_substeps )  THEN
518             CALL pres
519          ENDIF
520
521!
522!--       If required, compute liquid water content
523          IF ( cloud_physics )  THEN
524             CALL calc_liquid_water_content
525             !$acc update device( ql )
526          ENDIF
527!
528!--       If required, compute virtual potential temperature
529          IF ( humidity )  THEN
530             CALL compute_vpt 
531             !$acc update device( vpt )
532          ENDIF 
533!
534!--       Compute the diffusion quantities
535          IF ( .NOT. constant_diffusion )  THEN
536
537!
538!--          First the vertical fluxes in the Prandtl layer are being computed
539             IF ( prandtl_layer )  THEN
540                CALL cpu_log( log_point(19), 'prandtl_fluxes', 'start' )
541                CALL prandtl_fluxes
542                CALL cpu_log( log_point(19), 'prandtl_fluxes', 'stop' )
543             ENDIF
544!
545!--          Compute the diffusion coefficients
546             CALL cpu_log( log_point(17), 'diffusivities', 'start' )
547             IF ( .NOT. humidity ) THEN
548                IF ( ocean )  THEN
549                   CALL diffusivities( prho, prho_reference )
550                ELSE
551                   CALL diffusivities( pt, pt_reference )
552                ENDIF
553             ELSE
554                CALL diffusivities( vpt, pt_reference )
555             ENDIF
556             CALL cpu_log( log_point(17), 'diffusivities', 'stop' )
557
558          ENDIF
559
560       ENDDO   ! Intermediate step loop
561
562!
563!--    Increase simulation time and output times
564       nr_timesteps_this_run      = nr_timesteps_this_run + 1
565       current_timestep_number    = current_timestep_number + 1
566       simulated_time             = simulated_time   + dt_3d
567       simulated_time_chr         = time_to_string( simulated_time )
568       time_since_reference_point = simulated_time - coupling_start_time
569
570       IF ( simulated_time >= skip_time_data_output_av )  THEN
571          time_do_av         = time_do_av       + dt_3d
572       ENDIF
573       IF ( simulated_time >= skip_time_do2d_xy )  THEN
574          time_do2d_xy       = time_do2d_xy     + dt_3d
575       ENDIF
576       IF ( simulated_time >= skip_time_do2d_xz )  THEN
577          time_do2d_xz       = time_do2d_xz     + dt_3d
578       ENDIF
579       IF ( simulated_time >= skip_time_do2d_yz )  THEN
580          time_do2d_yz       = time_do2d_yz     + dt_3d
581       ENDIF
582       IF ( simulated_time >= skip_time_do3d    )  THEN
583          time_do3d          = time_do3d        + dt_3d
584       ENDIF
585       DO  mid = 1, masks
586          IF ( simulated_time >= skip_time_domask(mid) )  THEN
587             time_domask(mid)= time_domask(mid) + dt_3d
588          ENDIF
589       ENDDO
590       time_dvrp          = time_dvrp        + dt_3d
591       IF ( simulated_time >= skip_time_dosp )  THEN
592          time_dosp       = time_dosp        + dt_3d
593       ENDIF
594       time_dots          = time_dots        + dt_3d
595       IF ( .NOT. first_call_lpm )  THEN
596          time_dopts      = time_dopts       + dt_3d
597       ENDIF
598       IF ( simulated_time >= skip_time_dopr )  THEN
599          time_dopr       = time_dopr        + dt_3d
600       ENDIF
601       time_dopr_listing          = time_dopr_listing        + dt_3d
602       time_run_control   = time_run_control + dt_3d
603
604!
605!--    Data exchange between coupled models
606       IF ( coupling_mode /= 'uncoupled'  .AND.  run_coupled )  THEN
607          time_coupling = time_coupling + dt_3d
608
609!
610!--       In case of model termination initiated by the local model
611!--       (terminate_coupled > 0), the coupler must be skipped because it would
612!--       cause an MPI intercomminucation hang.
613!--       If necessary, the coupler will be called at the beginning of the
614!--       next restart run.
615          DO WHILE ( time_coupling >= dt_coupling .AND. terminate_coupled == 0 )
616             CALL surface_coupler
617             time_coupling = time_coupling - dt_coupling
618          ENDDO
619       ENDIF
620
621!
622!--    Execute user-defined actions
623       CALL user_actions( 'after_integration' )
624
625!
626!--    If Galilei transformation is used, determine the distance that the
627!--    model has moved so far
628       IF ( galilei_transformation )  THEN
629          advected_distance_x = advected_distance_x + u_gtrans * dt_3d
630          advected_distance_y = advected_distance_y + v_gtrans * dt_3d
631       ENDIF
632
633!
634!--    Check, if restart is necessary (because cpu-time is expiring or
635!--    because it is forced by user) and set stop flag
636!--    This call is skipped if the remote model has already initiated a restart.
637       IF ( .NOT. terminate_run )  CALL check_for_restart
638
639!
640!--    Carry out statistical analysis and output at the requested output times.
641!--    The MOD function is used for calculating the output time counters (like
642!--    time_dopr) in order to regard a possible decrease of the output time
643!--    interval in case of restart runs
644
645!
646!--    Set a flag indicating that so far no statistics have been created
647!--    for this time step
648       flow_statistics_called = .FALSE.
649
650!
651!--    If required, call flow_statistics for averaging in time
652       IF ( averaging_interval_pr /= 0.0  .AND.  &
653            ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.  &
654            simulated_time >= skip_time_dopr )  THEN
655          time_dopr_av = time_dopr_av + dt_3d
656          IF ( time_dopr_av >= dt_averaging_input_pr )  THEN
657             do_sum = .TRUE.
658             time_dopr_av = MOD( time_dopr_av, &
659                                    MAX( dt_averaging_input_pr, dt_3d ) )
660          ENDIF
661       ENDIF
662       IF ( do_sum )  CALL flow_statistics
663
664!
665!--    Sum-up 3d-arrays for later output of time-averaged 2d/3d/masked data
666       IF ( averaging_interval /= 0.0  .AND.                                &
667            ( dt_data_output_av - time_do_av ) <= averaging_interval  .AND. &
668            simulated_time >= skip_time_data_output_av )                    &
669       THEN
670          time_do_sla = time_do_sla + dt_3d
671          IF ( time_do_sla >= dt_averaging_input )  THEN
672             CALL sum_up_3d_data
673             average_count_3d = average_count_3d + 1
674             time_do_sla = MOD( time_do_sla, MAX( dt_averaging_input, dt_3d ) )
675          ENDIF
676       ENDIF
677
678!
679!--    Calculate spectra for time averaging
680       IF ( averaging_interval_sp /= 0.0  .AND.  &
681            ( dt_dosp - time_dosp ) <= averaging_interval_sp  .AND.  &
682            simulated_time >= skip_time_dosp )  THEN
683          time_dosp_av = time_dosp_av + dt_3d
684          IF ( time_dosp_av >= dt_averaging_input_pr )  THEN
685             CALL calc_spectra
686             time_dosp_av = MOD( time_dosp_av, &
687                                 MAX( dt_averaging_input_pr, dt_3d ) )
688          ENDIF
689       ENDIF
690
691!
692!--    Computation and output of run control parameters.
693!--    This is also done whenever perturbations have been imposed
694       IF ( time_run_control >= dt_run_control  .OR.                     &
695            timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created ) &
696       THEN
697          CALL run_control
698          IF ( time_run_control >= dt_run_control )  THEN
699             time_run_control = MOD( time_run_control, &
700                                     MAX( dt_run_control, dt_3d ) )
701          ENDIF
702       ENDIF
703
704!
705!--    Profile output (ASCII) on file
706       IF ( time_dopr_listing >= dt_dopr_listing )  THEN
707          CALL print_1d
708          time_dopr_listing = MOD( time_dopr_listing, MAX( dt_dopr_listing, &
709                                                           dt_3d ) )
710       ENDIF
711
712!
713!--    Graphic output for PROFIL
714       IF ( time_dopr >= dt_dopr )  THEN
715          IF ( dopr_n /= 0 )  CALL data_output_profiles
716          time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
717          time_dopr_av = 0.0    ! due to averaging (see above)
718       ENDIF
719
720!
721!--    Graphic output for time series
722       IF ( time_dots >= dt_dots )  THEN
723          CALL data_output_tseries
724          time_dots = MOD( time_dots, MAX( dt_dots, dt_3d ) )
725       ENDIF
726
727!
728!--    Output of spectra (formatted for use with PROFIL), in case of no
729!--    time averaging, spectra has to be calculated before
730       IF ( time_dosp >= dt_dosp )  THEN
731          IF ( average_count_sp == 0 )  CALL calc_spectra
732          CALL data_output_spectra
733          time_dosp = MOD( time_dosp, MAX( dt_dosp, dt_3d ) )
734       ENDIF
735
736!
737!--    2d-data output (cross-sections)
738       IF ( time_do2d_xy >= dt_do2d_xy )  THEN
739          CALL data_output_2d( 'xy', 0 )
740          time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
741       ENDIF
742       IF ( time_do2d_xz >= dt_do2d_xz )  THEN
743          CALL data_output_2d( 'xz', 0 )
744          time_do2d_xz = MOD( time_do2d_xz, MAX( dt_do2d_xz, dt_3d ) )
745       ENDIF
746       IF ( time_do2d_yz >= dt_do2d_yz )  THEN
747          CALL data_output_2d( 'yz', 0 )
748          time_do2d_yz = MOD( time_do2d_yz, MAX( dt_do2d_yz, dt_3d ) )
749       ENDIF
750
751!
752!--    3d-data output (volume data)
753       IF ( time_do3d >= dt_do3d )  THEN
754          CALL data_output_3d( 0 )
755          time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
756       ENDIF
757
758!
759!--    masked data output
760       DO  mid = 1, masks
761          IF ( time_domask(mid) >= dt_domask(mid) )  THEN
762             CALL data_output_mask( 0 )
763             time_domask(mid) = MOD( time_domask(mid),  &
764                                     MAX( dt_domask(mid), dt_3d ) )
765          ENDIF
766       ENDDO
767
768!
769!--    Output of time-averaged 2d/3d/masked data
770       IF ( time_do_av >= dt_data_output_av )  THEN
771          CALL average_3d_data
772          CALL data_output_2d( 'xy', 1 )
773          CALL data_output_2d( 'xz', 1 )
774          CALL data_output_2d( 'yz', 1 )
775          CALL data_output_3d( 1 )
776          DO  mid = 1, masks
777             CALL data_output_mask( 1 )
778          ENDDO
779          time_do_av = MOD( time_do_av, MAX( dt_data_output_av, dt_3d ) )
780       ENDIF
781
782!
783!--    Output of particle time series
784       IF ( particle_advection )  THEN
785          IF ( time_dopts >= dt_dopts  .OR. &
786               ( simulated_time >= particle_advection_start  .AND. &
787                 first_call_lpm ) )  THEN
788             CALL data_output_ptseries
789             time_dopts = MOD( time_dopts, MAX( dt_dopts, dt_3d ) )
790          ENDIF
791       ENDIF
792
793!
794!--    Output of dvrp-graphics (isosurface, particles, slicer)
795#if defined( __dvrp_graphics )
796       CALL DVRP_LOG_EVENT( -2, current_timestep_number-1 )
797#endif
798       IF ( time_dvrp >= dt_dvrp )  THEN
799          CALL data_output_dvrp
800          time_dvrp = MOD( time_dvrp, MAX( dt_dvrp, dt_3d ) )
801       ENDIF
802#if defined( __dvrp_graphics )
803       CALL DVRP_LOG_EVENT( 2, current_timestep_number )
804#endif
805
806!
807!--    If required, set the heat flux for the next time step at a random value
808       IF ( constant_heatflux  .AND.  random_heatflux )  CALL disturb_heatflux
809
810!
811!--    Execute user-defined actions
812       CALL user_actions( 'after_timestep' )
813
814       CALL cpu_log( log_point_s(10), 'timesteps', 'stop' )
815
816
817    ENDDO   ! time loop
818
819#if defined( __dvrp_graphics )
820    CALL DVRP_LOG_EVENT( -2, current_timestep_number )
821#endif
822
823 END SUBROUTINE time_integration
Note: See TracBrowser for help on using the repository browser.