Home

Context Navigation

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

Last change on this file since 1300 was 1277, checked in by heinze, 11 years ago
Last commit documented
Property svn:keywords set to `Id`
File size: 32.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |