Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |