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

Last change on this file since 2232 was 2232, checked in by suehring, 4 years ago

Adjustments according new topography and surface-modelling concept implemented

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