Home

Context Navigation

time_integration.f90 @ 1180

Last change on this file since 1180 was 1179, checked in by raasch, 12 years ago

New:
---
Initial profiles can be used as reference state in the buoyancy term. New parameter
reference_state introduced. Calculation and handling of reference state in buoyancy term revised.
binary version for restart files changed from 3.9 to 3.9a (no downward compatibility!),
initial profile for rho added to hom (id=77)

Errors:

small bugfix for background communication (time_integration)

Property svn:keywords set to Id

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

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

Context Navigation

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

Download in other formats: