Home

Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |