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

Last change on this file since 1015 was 1015, checked in by raasch, 9 years ago

Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added.
Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed.
mask array is set zero for ghost boundaries

  • Property svn:keywords set to Id
File size: 21.1 KB
Line 
1 SUBROUTINE time_integration
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6! +call of prognostic_equations_acc
7!
8! Former revisions:
9! -----------------
10! $Id: time_integration.f90 1015 2012-09-27 09:23:24Z raasch $
11!
12! 1001 2012-09-13 14:08:46Z raasch
13! all actions concerning leapfrog- and upstream-spline-scheme removed
14!
15! 849 2012-03-15 10:35:09Z raasch
16! advec_particles renamed lpm, first_call_advec_particles renamed first_call_lpm
17!
18! 825 2012-02-19 03:03:44Z raasch
19! wang_collision_kernel renamed wang_kernel
20!
21! 790 2011-11-29 03:11:20Z raasch
22! exchange of ghostpoints for array diss
23!
24! 707 2011-03-29 11:39:40Z raasch
25! bc_lr/ns replaced by bc_lr/ns_cyc, calls of exchange_horiz are modified,
26! adaption to sloping surface
27!
28! 667  2010-12-23 12:06:00Z suehring/gryschka
29! Calls of exchange_horiz are modified.
30! Adaption to slooping surface.
31!
32! 449 2010-02-02 11:23:59Z raasch
33! Bugfix: exchange of ghost points for prho included
34!
35! 410 2009-12-04 17:05:40Z letzel
36! masked data output
37!
38! 388 2009-09-23 09:40:33Z raasch
39! Using prho instead of rho in diffusvities.
40! Coupling with independent precursor runs.
41! Bugfix: output of particle time series only if particle advection is switched
42!         on
43!
44! 151 2008-03-07 13:42:18Z raasch
45! inflow turbulence is imposed by calling new routine inflow_turbulence
46!
47! 108 2007-08-24 15:10:38Z letzel
48! Call of new routine surface_coupler,
49! presure solver is called after the first Runge-Kutta substep instead of the
50! last in case that call_psolver_at_all_substeps = .F.; for this case, the
51! random perturbation has to be added to the velocity fields also after the
52! first substep
53!
54! 97 2007-06-21 08:23:15Z raasch
55! diffusivities is called with argument rho in case of ocean runs,
56! new argument pt_/prho_reference in calls of diffusivities,
57! ghostpoint exchange for salinity and density
58!
59! 87 2007-05-22 15:46:47Z raasch
60! var_hom renamed pr_palm
61!
62! 75 2007-03-22 09:54:05Z raasch
63! Move call of user_actions( 'after_integration' ) below increment of times
64! and counters,
65! calls of prognostic_equations_.. changed to .._noopt, .._cache, and
66! .._vector, these calls are now controlled by switch loop_optimization,
67! uxrp, vynp eliminated, 2nd+3rd argument removed from exchange horiz,
68! moisture renamed humidity
69!
70! RCS Log replace by Id keyword, revision history cleaned up
71!
72! Revision 1.8  2006/08/22 14:16:05  raasch
73! Disturbances are imposed only for the last Runge-Kutta-substep
74!
75! Revision 1.2  2004/04/30 13:03:40  raasch
76! decalpha-specific warning removed, routine name changed to time_integration,
77! particle advection is carried out only once during the intermediate steps,
78! impulse_advec renamed momentum_advec
79!
80! Revision 1.1  1997/08/11 06:19:04  raasch
81! Initial revision
82!
83!
84! Description:
85! ------------
86! Integration in time of the model equations, statistical analysis and graphic
87! output
88!------------------------------------------------------------------------------!
89
90    USE arrays_3d
91    USE averaging
92    USE control_parameters
93    USE cpulog
94#if defined( __dvrp_graphics )
95    USE DVRP
96#endif
97    USE grid_variables
98    USE indices
99    USE interaction_droplets_ptq_mod
100    USE interfaces
101    USE particle_attributes
102    USE pegrid
103    USE prognostic_equations_mod
104    USE statistics
105    USE user_actions_mod
106
107    IMPLICIT NONE
108
109    CHARACTER (LEN=9) ::  time_to_string
110    INTEGER ::  i, j, k
111
112!
113!-- At the beginning of a simulation determine the time step as well as
114!-- determine and print out the run control parameters
115    IF ( simulated_time == 0.0 )  CALL timestep
116
117    CALL run_control
118
119
120!
121!-- Data exchange between coupled models in case that a call has been omitted
122!-- at the end of the previous run of a job chain.
123    IF ( coupling_mode /= 'uncoupled'  .AND.  run_coupled )  THEN
124!
125!--    In case of model termination initiated by the local model the coupler
126!--    must not be called because this would again cause an MPI hang.
127       DO WHILE ( time_coupling >= dt_coupling .AND. terminate_coupled == 0 )
128          CALL surface_coupler
129          time_coupling = time_coupling - dt_coupling
130       ENDDO
131       IF (time_coupling == 0.0 .AND. time_since_reference_point < dt_coupling)&
132       THEN
133          time_coupling = time_since_reference_point
134       ENDIF
135    ENDIF
136
137
138#if defined( __dvrp_graphics )
139!
140!-- Time measurement with dvrp software 
141    CALL DVRP_LOG_EVENT( 2, current_timestep_number )
142#endif
143
144!
145!-- Start of the time loop
146    DO  WHILE ( simulated_time < end_time  .AND.  .NOT. stop_dt  .AND. &
147                .NOT. terminate_run )
148
149       CALL cpu_log( log_point_s(10), 'timesteps', 'start' )
150!
151!--    Determine size of next time step
152       IF ( simulated_time /= 0.0 )  CALL timestep
153!
154!--    Execute the user-defined actions
155       CALL user_actions( 'before_timestep' )
156
157!
158!--    Start of intermediate step loop
159       intermediate_timestep_count = 0
160       DO  WHILE ( intermediate_timestep_count < &
161                   intermediate_timestep_count_max )
162
163          intermediate_timestep_count = intermediate_timestep_count + 1
164
165!
166!--       Set the steering factors for the prognostic equations which depend
167!--       on the timestep scheme
168          CALL timestep_scheme_steering
169
170!
171!--       Solve the prognostic equations. A fast cache optimized version with
172!--       only one single loop is used in case of Piascek-Williams advection
173!--       scheme. NEC vector machines use a different version, because
174!--       in the other versions a good vectorization is prohibited due to
175!--       inlining problems.
176          IF ( loop_optimization == 'vector' )  THEN
177             CALL prognostic_equations_vector
178          ELSEIF ( loop_optimization == 'acc' )  THEN
179             CALL prognostic_equations_acc
180          ELSE
181             IF ( scalar_advec == 'bc-scheme' )  THEN
182                CALL prognostic_equations_noopt
183             ELSE
184                CALL prognostic_equations_cache
185             ENDIF
186          ENDIF
187
188!
189!--       Particle transport/physics with the Lagrangian particle model
190!--       (only once during intermediate steps, because it uses an Euler-step)
191          IF ( particle_advection  .AND.                         &
192               simulated_time >= particle_advection_start  .AND. &
193               intermediate_timestep_count == 1 )  THEN
194             CALL lpm
195             first_call_lpm = .FALSE.
196          ENDIF
197
198!
199!--       Interaction of droplets with temperature and specific humidity.
200!--       Droplet condensation and evaporation is calculated within
201!--       advec_particles.
202          IF ( cloud_droplets  .AND.  &
203               intermediate_timestep_count == intermediate_timestep_count_max )&
204          THEN
205             CALL interaction_droplets_ptq
206          ENDIF
207
208!
209!--       Exchange of ghost points (lateral boundary conditions)
210          CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'start' )
211          CALL exchange_horiz( u_p, nbgp )
212          CALL exchange_horiz( v_p, nbgp )
213          CALL exchange_horiz( w_p, nbgp )
214          CALL exchange_horiz( pt_p, nbgp )
215          IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e_p, nbgp )
216          IF ( ocean )  THEN
217             CALL exchange_horiz( sa_p, nbgp )
218             CALL exchange_horiz( rho, nbgp )
219             CALL exchange_horiz( prho, nbgp )
220          ENDIF
221          IF (humidity  .OR.  passive_scalar)  CALL exchange_horiz( q_p, nbgp )
222          IF ( cloud_droplets )  THEN
223             CALL exchange_horiz( ql, nbgp )
224             CALL exchange_horiz( ql_c, nbgp )
225             CALL exchange_horiz( ql_v, nbgp )
226             CALL exchange_horiz( ql_vp, nbgp )
227          ENDIF
228          IF ( wang_kernel )  CALL exchange_horiz( diss, nbgp )
229
230          CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
231
232!
233!--       Boundary conditions for the prognostic quantities (except of the
234!--       velocities at the outflow in case of a non-cyclic lateral wall)
235          CALL boundary_conds( 'main' )
236
237!
238!--       Swap the time levels in preparation for the next time step.
239          CALL swap_timelevel
240
241!
242!--       Temperature offset must be imposed at cyclic boundaries in x-direction
243!--       when a sloping surface is used
244          IF ( sloping_surface )  THEN
245             IF ( nxl ==  0 )  pt(:,:,nxlg:nxl-1) = pt(:,:,nxlg:nxl-1) - &
246                                                    pt_slope_offset
247             IF ( nxr == nx )  pt(:,:,nxr+1:nxrg) = pt(:,:,nxr+1:nxrg) + &
248                                                    pt_slope_offset
249          ENDIF
250
251!
252!--       Impose a turbulent inflow using the recycling method
253          IF ( turbulent_inflow )  CALL  inflow_turbulence
254
255!
256!--       Impose a random perturbation on the horizontal velocity field
257          IF ( create_disturbances  .AND.                                      &
258               ( call_psolver_at_all_substeps  .AND.                           &
259               intermediate_timestep_count == intermediate_timestep_count_max )&
260          .OR. ( .NOT. call_psolver_at_all_substeps  .AND.                     &
261               intermediate_timestep_count == 1 ) )                            &
262          THEN
263             time_disturb = time_disturb + dt_3d
264             IF ( time_disturb >= dt_disturb )  THEN
265                IF ( hom(nzb+5,1,pr_palm,0) < disturbance_energy_limit )  THEN
266                   CALL disturb_field( nzb_u_inner, tend, u )
267                   CALL disturb_field( nzb_v_inner, tend, v )
268                ELSEIF ( .NOT. bc_lr_cyc  .OR.  .NOT. bc_ns_cyc )  THEN
269!
270!--                Runs with a non-cyclic lateral wall need perturbations
271!--                near the inflow throughout the whole simulation
272                   dist_range = 1
273                   CALL disturb_field( nzb_u_inner, tend, u )
274                   CALL disturb_field( nzb_v_inner, tend, v )
275                   dist_range = 0
276                ENDIF
277                time_disturb = time_disturb - dt_disturb
278             ENDIF
279          ENDIF
280
281!
282!--       Reduce the velocity divergence via the equation for perturbation
283!--       pressure.
284          IF ( intermediate_timestep_count == 1  .OR. &
285                call_psolver_at_all_substeps )  THEN
286             CALL pres
287          ENDIF
288!
289!--       Update device memory for calculating diffusion quantities and for next
290!--       timestep
291          !$acc update device( e, pt, u, v, w )
292          !$acc update device( q )  if ( allocated( q ) )
293
294!
295!--       If required, compute virtuell potential temperature
296          IF ( humidity )  THEN
297             CALL compute_vpt
298             !$acc update device( vpt )
299          ENDIF
300
301!
302!--       If required, compute liquid water content
303          IF ( cloud_physics )  THEN
304             CALL calc_liquid_water_content
305             !$acc update device( ql )
306          ENDIF
307
308!
309!--       Compute the diffusion quantities
310          IF ( .NOT. constant_diffusion )  THEN
311
312!
313!--          First the vertical fluxes in the Prandtl layer are being computed
314             IF ( prandtl_layer )  THEN
315                CALL cpu_log( log_point(19), 'prandtl_fluxes', 'start' )
316                CALL prandtl_fluxes
317                CALL cpu_log( log_point(19), 'prandtl_fluxes', 'stop' )
318!
319!++             Statistics still require updates on host
320                !$acc update host( qs, qsws, rif, shf, ts )
321             ENDIF
322
323!
324!--          Compute the diffusion coefficients
325             CALL cpu_log( log_point(17), 'diffusivities', 'start' )
326             IF ( .NOT. humidity ) THEN
327                IF ( ocean )  THEN
328                   CALL diffusivities( prho, prho_reference )
329                ELSE
330                   CALL diffusivities( pt, pt_reference )
331                ENDIF
332             ELSE
333                CALL diffusivities( vpt, pt_reference )
334             ENDIF
335             CALL cpu_log( log_point(17), 'diffusivities', 'stop' )
336!
337!++          Statistics still require update of diffusivities on host
338             !$acc update host( kh, km )
339
340          ENDIF
341
342       ENDDO   ! Intermediate step loop
343
344!
345!--    Increase simulation time and output times
346       current_timestep_number    = current_timestep_number + 1
347       simulated_time             = simulated_time   + dt_3d
348       simulated_time_chr         = time_to_string( simulated_time )
349       time_since_reference_point = simulated_time - coupling_start_time
350
351       IF ( simulated_time >= skip_time_data_output_av )  THEN
352          time_do_av         = time_do_av       + dt_3d
353       ENDIF
354       IF ( simulated_time >= skip_time_do2d_xy )  THEN
355          time_do2d_xy       = time_do2d_xy     + dt_3d
356       ENDIF
357       IF ( simulated_time >= skip_time_do2d_xz )  THEN
358          time_do2d_xz       = time_do2d_xz     + dt_3d
359       ENDIF
360       IF ( simulated_time >= skip_time_do2d_yz )  THEN
361          time_do2d_yz       = time_do2d_yz     + dt_3d
362       ENDIF
363       IF ( simulated_time >= skip_time_do3d    )  THEN
364          time_do3d          = time_do3d        + dt_3d
365       ENDIF
366       DO  mid = 1, masks
367          IF ( simulated_time >= skip_time_domask(mid) )  THEN
368             time_domask(mid)= time_domask(mid) + dt_3d
369          ENDIF
370       ENDDO
371       time_dvrp          = time_dvrp        + dt_3d
372       IF ( simulated_time >= skip_time_dosp )  THEN
373          time_dosp       = time_dosp        + dt_3d
374       ENDIF
375       time_dots          = time_dots        + dt_3d
376       IF ( .NOT. first_call_lpm )  THEN
377          time_dopts      = time_dopts       + dt_3d
378       ENDIF
379       IF ( simulated_time >= skip_time_dopr )  THEN
380          time_dopr       = time_dopr        + dt_3d
381       ENDIF
382       time_dopr_listing          = time_dopr_listing        + dt_3d
383       time_run_control   = time_run_control + dt_3d
384
385!
386!--    Data exchange between coupled models
387       IF ( coupling_mode /= 'uncoupled'  .AND.  run_coupled )  THEN
388          time_coupling = time_coupling + dt_3d
389
390!
391!--       In case of model termination initiated by the local model
392!--       (terminate_coupled > 0), the coupler must be skipped because it would
393!--       cause an MPI intercomminucation hang.
394!--       If necessary, the coupler will be called at the beginning of the
395!--       next restart run.
396          DO WHILE ( time_coupling >= dt_coupling .AND. terminate_coupled == 0 )
397             CALL surface_coupler
398             time_coupling = time_coupling - dt_coupling
399          ENDDO
400       ENDIF
401
402!
403!--    Execute user-defined actions
404       CALL user_actions( 'after_integration' )
405
406!
407!--    If Galilei transformation is used, determine the distance that the
408!--    model has moved so far
409       IF ( galilei_transformation )  THEN
410          advected_distance_x = advected_distance_x + u_gtrans * dt_3d
411          advected_distance_y = advected_distance_y + v_gtrans * dt_3d
412       ENDIF
413
414!
415!--    Check, if restart is necessary (because cpu-time is expiring or
416!--    because it is forced by user) and set stop flag
417!--    This call is skipped if the remote model has already initiated a restart.
418       IF ( .NOT. terminate_run )  CALL check_for_restart
419
420!
421!--    Carry out statistical analysis and output at the requested output times.
422!--    The MOD function is used for calculating the output time counters (like
423!--    time_dopr) in order to regard a possible decrease of the output time
424!--    interval in case of restart runs
425
426!
427!--    Set a flag indicating that so far no statistics have been created
428!--    for this time step
429       flow_statistics_called = .FALSE.
430
431!
432!--    If required, call flow_statistics for averaging in time
433       IF ( averaging_interval_pr /= 0.0  .AND.  &
434            ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.  &
435            simulated_time >= skip_time_dopr )  THEN
436          time_dopr_av = time_dopr_av + dt_3d
437          IF ( time_dopr_av >= dt_averaging_input_pr )  THEN
438             do_sum = .TRUE.
439             time_dopr_av = MOD( time_dopr_av, &
440                                    MAX( dt_averaging_input_pr, dt_3d ) )
441          ENDIF
442       ENDIF
443       IF ( do_sum )  CALL flow_statistics
444
445!
446!--    Sum-up 3d-arrays for later output of time-averaged 2d/3d/masked data
447       IF ( averaging_interval /= 0.0  .AND.                                &
448            ( dt_data_output_av - time_do_av ) <= averaging_interval  .AND. &
449            simulated_time >= skip_time_data_output_av )                    &
450       THEN
451          time_do_sla = time_do_sla + dt_3d
452          IF ( time_do_sla >= dt_averaging_input )  THEN
453             CALL sum_up_3d_data
454             average_count_3d = average_count_3d + 1
455             time_do_sla = MOD( time_do_sla, MAX( dt_averaging_input, dt_3d ) )
456          ENDIF
457       ENDIF
458
459!
460!--    Calculate spectra for time averaging
461       IF ( averaging_interval_sp /= 0.0  .AND.  &
462            ( dt_dosp - time_dosp ) <= averaging_interval_sp  .AND.  &
463            simulated_time >= skip_time_dosp )  THEN
464          time_dosp_av = time_dosp_av + dt_3d
465          IF ( time_dosp_av >= dt_averaging_input_pr )  THEN
466             CALL calc_spectra
467             time_dosp_av = MOD( time_dosp_av, &
468                                 MAX( dt_averaging_input_pr, dt_3d ) )
469          ENDIF
470       ENDIF
471
472!
473!--    Computation and output of run control parameters.
474!--    This is also done whenever perturbations have been imposed
475       IF ( time_run_control >= dt_run_control  .OR.                     &
476            timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created ) &
477       THEN
478          CALL run_control
479          IF ( time_run_control >= dt_run_control )  THEN
480             time_run_control = MOD( time_run_control, &
481                                     MAX( dt_run_control, dt_3d ) )
482          ENDIF
483       ENDIF
484
485!
486!--    Profile output (ASCII) on file
487       IF ( time_dopr_listing >= dt_dopr_listing )  THEN
488          CALL print_1d
489          time_dopr_listing = MOD( time_dopr_listing, MAX( dt_dopr_listing, &
490                                                           dt_3d ) )
491       ENDIF
492
493!
494!--    Graphic output for PROFIL
495       IF ( time_dopr >= dt_dopr )  THEN
496          IF ( dopr_n /= 0 )  CALL data_output_profiles
497          time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
498          time_dopr_av = 0.0    ! due to averaging (see above)
499       ENDIF
500
501!
502!--    Graphic output for time series
503       IF ( time_dots >= dt_dots )  THEN
504          CALL data_output_tseries
505          time_dots = MOD( time_dots, MAX( dt_dots, dt_3d ) )
506       ENDIF
507
508!
509!--    Output of spectra (formatted for use with PROFIL), in case of no
510!--    time averaging, spectra has to be calculated before
511       IF ( time_dosp >= dt_dosp )  THEN
512          IF ( average_count_sp == 0 )  CALL calc_spectra
513          CALL data_output_spectra
514          time_dosp = MOD( time_dosp, MAX( dt_dosp, dt_3d ) )
515       ENDIF
516
517!
518!--    2d-data output (cross-sections)
519       IF ( time_do2d_xy >= dt_do2d_xy )  THEN
520          CALL data_output_2d( 'xy', 0 )
521          time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
522       ENDIF
523       IF ( time_do2d_xz >= dt_do2d_xz )  THEN
524          CALL data_output_2d( 'xz', 0 )
525          time_do2d_xz = MOD( time_do2d_xz, MAX( dt_do2d_xz, dt_3d ) )
526       ENDIF
527       IF ( time_do2d_yz >= dt_do2d_yz )  THEN
528          CALL data_output_2d( 'yz', 0 )
529          time_do2d_yz = MOD( time_do2d_yz, MAX( dt_do2d_yz, dt_3d ) )
530       ENDIF
531
532!
533!--    3d-data output (volume data)
534       IF ( time_do3d >= dt_do3d )  THEN
535          CALL data_output_3d( 0 )
536          time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
537       ENDIF
538
539!
540!--    masked data output
541       DO  mid = 1, masks
542          IF ( time_domask(mid) >= dt_domask(mid) )  THEN
543             CALL data_output_mask( 0 )
544             time_domask(mid) = MOD( time_domask(mid),  &
545                                     MAX( dt_domask(mid), dt_3d ) )
546          ENDIF
547       ENDDO
548
549!
550!--    Output of time-averaged 2d/3d/masked data
551       IF ( time_do_av >= dt_data_output_av )  THEN
552          CALL average_3d_data
553          CALL data_output_2d( 'xy', 1 )
554          CALL data_output_2d( 'xz', 1 )
555          CALL data_output_2d( 'yz', 1 )
556          CALL data_output_3d( 1 )
557          DO  mid = 1, masks
558             CALL data_output_mask( 1 )
559          ENDDO
560          time_do_av = MOD( time_do_av, MAX( dt_data_output_av, dt_3d ) )
561       ENDIF
562
563!
564!--    Output of particle time series
565       IF ( particle_advection )  THEN
566          IF ( time_dopts >= dt_dopts  .OR. &
567               ( simulated_time >= particle_advection_start  .AND. &
568                 first_call_lpm ) )  THEN
569             CALL data_output_ptseries
570             time_dopts = MOD( time_dopts, MAX( dt_dopts, dt_3d ) )
571          ENDIF
572       ENDIF
573
574!
575!--    Output of dvrp-graphics (isosurface, particles, slicer)
576#if defined( __dvrp_graphics )
577       CALL DVRP_LOG_EVENT( -2, current_timestep_number-1 )
578#endif
579       IF ( time_dvrp >= dt_dvrp )  THEN
580          CALL data_output_dvrp
581          time_dvrp = MOD( time_dvrp, MAX( dt_dvrp, dt_3d ) )
582       ENDIF
583#if defined( __dvrp_graphics )
584       CALL DVRP_LOG_EVENT( 2, current_timestep_number )
585#endif
586
587!
588!--    If required, set the heat flux for the next time step at a random value
589       IF ( constant_heatflux  .AND.  random_heatflux )  CALL disturb_heatflux
590
591!
592!--    Execute user-defined actions
593       CALL user_actions( 'after_timestep' )
594
595       CALL cpu_log( log_point_s(10), 'timesteps', 'stop' )
596
597
598    ENDDO   ! time loop
599
600#if defined( __dvrp_graphics )
601    CALL DVRP_LOG_EVENT( -2, current_timestep_number )
602#endif
603
604 END SUBROUTINE time_integration
Note: See TracBrowser for help on using the repository browser.