Home

Context Navigation

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

Last change on this file since 1463 was 1403, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 37.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |