Home

Context Navigation

source: palm/trunk/SOURCE/prognostic_equations.f90 @ 1371

Last change on this file since 1371 was 1366, checked in by boeske, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 82.7 KB

Rev	Line
[736]	1	MODULE prognostic_equations_mod
	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	!
[1310]	17	! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]	18	!--------------------------------------------------------------------------------!
	19	!
[736]	20	! Current revisions:
[1092]	21	! ------------------
[1366]	22	!
	23	!
	24	! Former revisions:
	25	! -----------------
	26	! $Id: prognostic_equations.f90 1366 2014-04-22 15:06:33Z raasch $
	27	!
	28	! 1365 2014-04-22 15:03:56Z boeske
[1365]	29	! Calls of ls_advec for large scale advection added,
	30	! subroutine subsidence is only called if use_subsidence_tendencies = .F.,
	31	! new argument ls_index added to the calls of subsidence
	32	! +ls_index
	33	!
	34	! 1361 2014-04-16 15:17:48Z hoffmann
[1361]	35	! Two-moment microphysics moved to the start of prognostic equations. This makes
	36	! the 3d arrays for tend_q, tend_qr, tend_pt and tend_pt redundant.
	37	! Additionally, it is allowed to call the microphysics just once during the time
	38	! step (not at each sub-time step).
	39	!
	40	! Two-moment cloud physics added for vector and accelerator optimization.
	41	!
	42	! 1353 2014-04-08 15:21:23Z heinze
[1354]	43	! REAL constants provided with KIND-attribute
	44	!
[1353]	45	! 1337 2014-03-25 15:11:48Z heinze
	46	! Bugfix: REAL constants provided with KIND-attribute
	47	!
	48	! 1332 2014-03-25 11:59:43Z suehring
[1333]	49	! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke
	50	!
[1332]	51	! 1330 2014-03-24 17:29:32Z suehring
[1331]	52	! In case of SGS-particle velocity advection of TKE is also allowed with
	53	! dissipative 5th-order scheme.
	54	!
[1321]	55	! 1320 2014-03-20 08:40:49Z raasch
[1320]	56	! ONLY-attribute added to USE-statements,
	57	! kind-parameters added to all INTEGER and REAL declaration statements,
	58	! kinds are defined in new module kinds,
	59	! old module precision_kind is removed,
	60	! revision history before 2012 removed,
	61	! comment fields (!:) to be used for variable explanations added to
	62	! all variable declaration statements
[1054]	63	!
[1319]	64	! 1318 2014-03-17 13:35:16Z raasch
	65	! module interfaces removed
	66	!
[1258]	67	! 1257 2013-11-08 15:18:40Z raasch
	68	! openacc loop vector clauses removed, independent clauses added
	69	!
[1247]	70	! 1246 2013-11-01 08:59:45Z heinze
	71	! enable nudging also for accelerator version
	72	!
[1242]	73	! 1241 2013-10-30 11:36:58Z heinze
	74	! usage of nudging enabled (so far not implemented for accelerator version)
	75	!
[1182]	76	! 1179 2013-06-14 05:57:58Z raasch
	77	! two arguments removed from routine buoyancy, ref_state updated on device
	78	!
[1132]	79	! 1128 2013-04-12 06:19:32Z raasch
	80	! those parts requiring global communication moved to time_integration,
	81	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
	82	! j_north
	83	!
[1116]	84	! 1115 2013-03-26 18:16:16Z hoffmann
	85	! optimized cloud physics: calculation of microphysical tendencies transfered
	86	! to microphysics.f90; qr and nr are only calculated if precipitation is required
	87	!
[1112]	88	! 1111 2013-03-08 23:54:10Z raasch
	89	! update directives for prognostic quantities removed
	90	!
[1107]	91	! 1106 2013-03-04 05:31:38Z raasch
	92	! small changes in code formatting
	93	!
[1093]	94	! 1092 2013-02-02 11:24:22Z raasch
	95	! unused variables removed
	96	!
[1054]	97	! 1053 2012-11-13 17:11:03Z hoffmann
[1053]	98	! implementation of two new prognostic equations for rain drop concentration (nr)
	99	! and rain water content (qr)
[979]	100	!
[1053]	101	! currently, only available for cache loop optimization
[1020]	102	!
[1037]	103	! 1036 2012-10-22 13:43:42Z raasch
	104	! code put under GPL (PALM 3.9)
	105	!
[1020]	106	! 1019 2012-09-28 06:46:45Z raasch
	107	! non-optimized version of prognostic_equations removed
	108	!
[1017]	109	! 1015 2012-09-27 09:23:24Z raasch
	110	! new branch prognostic_equations_acc
	111	! OpenACC statements added + code changes required for GPU optimization
	112	!
[1002]	113	! 1001 2012-09-13 14:08:46Z raasch
	114	! all actions concerning leapfrog- and upstream-spline-scheme removed
	115	!
[979]	116	! 978 2012-08-09 08:28:32Z fricke
[978]	117	! km_damp_x and km_damp_y removed in calls of diffusion_u and diffusion_v
	118	! add ptdf_x, ptdf_y for damping the potential temperature at the inflow
	119	! boundary in case of non-cyclic lateral boundaries
	120	! Bugfix: first thread index changes for WS-scheme at the inflow
[736]	121	!
[941]	122	! 940 2012-07-09 14:31:00Z raasch
	123	! temperature equation can be switched off
	124	!
[736]	125	! Revision 1.1 2000/04/13 14:56:27 schroeter
	126	! Initial revision
	127	!
	128	!
	129	! Description:
	130	! ------------
	131	! Solving the prognostic equations.
	132	!------------------------------------------------------------------------------!
	133
[1320]	134	USE arrays_3d, &
	135	ONLY: diss_l_e, diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr, &
	136	diss_l_sa, diss_s_e, diss_s_nr, diss_s_pt, diss_s_q, &
	137	diss_s_qr, diss_s_sa, e, e_p, flux_s_e, flux_s_nr, flux_s_pt, &
	138	flux_s_q, flux_s_qr, flux_s_sa, flux_l_e, flux_l_nr, &
	139	flux_l_pt, flux_l_q, flux_l_qr, flux_l_sa, nr, nr_p, nrsws, &
	140	nrswst, pt, ptdf_x, ptdf_y, pt_init, pt_p, prho, q, q_init, &
	141	q_p, qsws, qswst, qr, qr_p, qrsws, qrswst, rdf, rdf_sc, rho, &
[1361]	142	sa, sa_init, sa_p, saswsb, saswst, shf, tend, te_m, tnr_m, &
	143	tpt_m, tq_m, tqr_m, tsa_m, tswst, tu_m, tv_m, tw_m, u, ug, u_p, &
	144	v, vg, vpt, v_p, w, w_p
[1320]	145
	146	USE control_parameters, &
[1361]	147	ONLY: call_microphysics_at_all_substeps, cloud_physics, &
	148	constant_diffusion, cthf, dp_external, &
[1320]	149	dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d, humidity, &
	150	icloud_scheme, inflow_l, intermediate_timestep_count, &
[1365]	151	intermediate_timestep_count_max, large_scale_forcing, &
	152	large_scale_subsidence, neutral, nudging, ocean, outflow_l, &
	153	outflow_s, passive_scalar, plant_canopy, precipitation, &
	154	prho_reference, prho_reference, prho_reference, pt_reference, &
	155	pt_reference, pt_reference, radiation, scalar_advec, &
	156	scalar_advec, simulated_time, sloping_surface, timestep_scheme, &
	157	tsc, use_subsidence_tendencies, use_upstream_for_tke, &
[1320]	158	use_upstream_for_tke, use_upstream_for_tke, wall_heatflux, &
	159	wall_nrflux, wall_qflux, wall_qflux, wall_qflux, wall_qrflux, &
	160	wall_salinityflux, ws_scheme_mom, ws_scheme_sca
[736]	161
[1320]	162	USE cpulog, &
	163	ONLY: cpu_log, log_point
[736]	164
[1320]	165	USE eqn_state_seawater_mod, &
	166	ONLY: eqn_state_seawater
	167
	168	USE indices, &
	169	ONLY: i_left, i_right, j_north, j_south, nxl, nxlu, nxr, nyn, nys, &
	170	nysv, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt
	171
	172	USE advec_ws, &
	173	ONLY: advec_s_ws, advec_s_ws_acc, advec_u_ws, advec_u_ws_acc, &
	174	advec_v_ws, advec_v_ws_acc, advec_w_ws, advec_w_ws_acc
	175
	176	USE advec_s_pw_mod, &
	177	ONLY: advec_s_pw
	178
	179	USE advec_s_up_mod, &
	180	ONLY: advec_s_up
	181
	182	USE advec_u_pw_mod, &
	183	ONLY: advec_u_pw
	184
	185	USE advec_u_up_mod, &
	186	ONLY: advec_u_up
	187
	188	USE advec_v_pw_mod, &
	189	ONLY: advec_v_pw
	190
	191	USE advec_v_up_mod, &
	192	ONLY: advec_v_up
	193
	194	USE advec_w_pw_mod, &
	195	ONLY: advec_w_pw
	196
	197	USE advec_w_up_mod, &
	198	ONLY: advec_w_up
	199
	200	USE buoyancy_mod, &
	201	ONLY: buoyancy, buoyancy_acc
	202
	203	USE calc_precipitation_mod, &
	204	ONLY: calc_precipitation
	205
	206	USE calc_radiation_mod, &
	207	ONLY: calc_radiation
	208
	209	USE coriolis_mod, &
	210	ONLY: coriolis, coriolis_acc
	211
	212	USE diffusion_e_mod, &
	213	ONLY: diffusion_e, diffusion_e_acc
	214
	215	USE diffusion_s_mod, &
	216	ONLY: diffusion_s, diffusion_s_acc
	217
	218	USE diffusion_u_mod, &
	219	ONLY: diffusion_u, diffusion_u_acc
	220
	221	USE diffusion_v_mod, &
	222	ONLY: diffusion_v, diffusion_v_acc
	223
	224	USE diffusion_w_mod, &
	225	ONLY: diffusion_w, diffusion_w_acc
	226
	227	USE impact_of_latent_heat_mod, &
	228	ONLY: impact_of_latent_heat
	229
	230	USE kinds
	231
[1365]	232	USE ls_forcing_mod, &
	233	ONLY: ls_advec
	234
[1320]	235	USE microphysics_mod, &
	236	ONLY: microphysics_control
	237
	238	USE nudge_mod, &
	239	ONLY: nudge
	240
	241	USE plant_canopy_model_mod, &
	242	ONLY: plant_canopy_model
	243
	244	USE production_e_mod, &
	245	ONLY: production_e, production_e_acc
	246
	247	USE subsidence_mod, &
	248	ONLY: subsidence
	249
	250	USE user_actions_mod, &
	251	ONLY: user_actions
	252
	253
[736]	254	PRIVATE
[1019]	255	PUBLIC prognostic_equations_cache, prognostic_equations_vector, &
	256	prognostic_equations_acc
[736]	257
	258	INTERFACE prognostic_equations_cache
	259	MODULE PROCEDURE prognostic_equations_cache
	260	END INTERFACE prognostic_equations_cache
	261
	262	INTERFACE prognostic_equations_vector
	263	MODULE PROCEDURE prognostic_equations_vector
	264	END INTERFACE prognostic_equations_vector
	265
[1015]	266	INTERFACE prognostic_equations_acc
	267	MODULE PROCEDURE prognostic_equations_acc
	268	END INTERFACE prognostic_equations_acc
[736]	269
[1015]	270
[736]	271	CONTAINS
	272
	273
	274	SUBROUTINE prognostic_equations_cache
	275
	276	!------------------------------------------------------------------------------!
	277	! Version with one optimized loop over all equations. It is only allowed to
	278	! be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
	279	!
	280	! Here the calls of most subroutines are embedded in two DO loops over i and j,
	281	! so communication between CPUs is not allowed (does not make sense) within
	282	! these loops.
	283	!
	284	! (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
	285	!------------------------------------------------------------------------------!
	286
	287	IMPLICIT NONE
	288
[1320]	289	INTEGER(iwp) :: i !:
	290	INTEGER(iwp) :: i_omp_start !:
	291	INTEGER(iwp) :: j !:
	292	INTEGER(iwp) :: k !:
	293	INTEGER(iwp) :: omp_get_thread_num !:
	294	INTEGER(iwp) :: tn = 0 !:
	295
	296	LOGICAL :: loop_start !:
[736]	297
	298
	299	!
	300	!-- Time measurement can only be performed for the whole set of equations
	301	CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
	302
	303	!
	304	!-- Loop over all prognostic equations
	305	!$OMP PARALLEL private (i,i_omp_start,j,k,loop_start,tn)
	306
	307	!$ tn = omp_get_thread_num()
	308	loop_start = .TRUE.
	309	!$OMP DO
	310	DO i = nxl, nxr
	311
	312	!
	313	!-- Store the first loop index. It differs for each thread and is required
	314	!-- later in advec_ws
	315	IF ( loop_start ) THEN
	316	loop_start = .FALSE.
	317	i_omp_start = i
	318	ENDIF
[1365]	319
[736]	320	DO j = nys, nyn
	321	!
[1361]	322	!-- If required, calculate cloud microphysical impacts (two-moment scheme)
	323	IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &
	324	( intermediate_timestep_count == 1 .OR. &
	325	call_microphysics_at_all_substeps ) &
	326	) THEN
	327	CALL microphysics_control( i, j )
	328	ENDIF
	329	!
[736]	330	!-- Tendency terms for u-velocity component
	331	IF ( .NOT. outflow_l .OR. i > nxl ) THEN
	332
[1337]	333	tend(:,j,i) = 0.0_wp
[1001]	334	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	335	IF ( ws_scheme_mom ) THEN
[978]	336	IF ( ( inflow_l .OR. outflow_l ) .AND. i_omp_start == nxl ) THEN
[736]	337	CALL advec_u_ws( i, j, i_omp_start + 1, tn )
	338	ELSE
	339	CALL advec_u_ws( i, j, i_omp_start, tn )
	340	ENDIF
	341	ELSE
	342	CALL advec_u_pw( i, j )
	343	ENDIF
	344	ELSE
	345	CALL advec_u_up( i, j )
	346	ENDIF
[1001]	347	CALL diffusion_u( i, j )
[736]	348	CALL coriolis( i, j, 1 )
[940]	349	IF ( sloping_surface .AND. .NOT. neutral ) THEN
[1179]	350	CALL buoyancy( i, j, pt, 1 )
[940]	351	ENDIF
[736]	352
	353	!
	354	!-- Drag by plant canopy
	355	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 1 )
	356
	357	!
	358	!-- External pressure gradient
	359	IF ( dp_external ) THEN
	360	DO k = dp_level_ind_b+1, nzt
	361	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
	362	ENDDO
	363	ENDIF
	364
[1241]	365	!
	366	!-- Nudging
	367	IF ( nudging ) CALL nudge( i, j, simulated_time, 'u' )
	368
[736]	369	CALL user_actions( i, j, 'u-tendency' )
	370	!
	371	!-- Prognostic equation for u-velocity component
	372	DO k = nzb_u_inner(j,i)+1, nzt
[1001]	373	u_p(k,j,i) = u(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	374	tsc(3) * tu_m(k,j,i) ) &
	375	- tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
[736]	376	ENDDO
	377
	378	!
	379	!-- Calculate tendencies for the next Runge-Kutta step
	380	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	381	IF ( intermediate_timestep_count == 1 ) THEN
	382	DO k = nzb_u_inner(j,i)+1, nzt
	383	tu_m(k,j,i) = tend(k,j,i)
	384	ENDDO
	385	ELSEIF ( intermediate_timestep_count < &
	386	intermediate_timestep_count_max ) THEN
	387	DO k = nzb_u_inner(j,i)+1, nzt
[1337]	388	tu_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tu_m(k,j,i)
[736]	389	ENDDO
	390	ENDIF
	391	ENDIF
	392
	393	ENDIF
	394
	395	!
	396	!-- Tendency terms for v-velocity component
	397	IF ( .NOT. outflow_s .OR. j > nys ) THEN
	398
[1337]	399	tend(:,j,i) = 0.0_wp
[1001]	400	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	401	IF ( ws_scheme_mom ) THEN
	402	CALL advec_v_ws( i, j, i_omp_start, tn )
	403	ELSE
	404	CALL advec_v_pw( i, j )
	405	ENDIF
	406	ELSE
	407	CALL advec_v_up( i, j )
	408	ENDIF
[1001]	409	CALL diffusion_v( i, j )
[736]	410	CALL coriolis( i, j, 2 )
	411
	412	!
	413	!-- Drag by plant canopy
	414	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 2 )
	415
	416	!
	417	!-- External pressure gradient
	418	IF ( dp_external ) THEN
	419	DO k = dp_level_ind_b+1, nzt
	420	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
	421	ENDDO
	422	ENDIF
	423
[1241]	424	!
	425	!-- Nudging
	426	IF ( nudging ) CALL nudge( i, j, simulated_time, 'v' )
	427
[736]	428	CALL user_actions( i, j, 'v-tendency' )
	429	!
	430	!-- Prognostic equation for v-velocity component
	431	DO k = nzb_v_inner(j,i)+1, nzt
[1001]	432	v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	433	tsc(3) * tv_m(k,j,i) ) &
	434	- tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
[736]	435	ENDDO
	436
	437	!
	438	!-- Calculate tendencies for the next Runge-Kutta step
	439	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	440	IF ( intermediate_timestep_count == 1 ) THEN
	441	DO k = nzb_v_inner(j,i)+1, nzt
	442	tv_m(k,j,i) = tend(k,j,i)
	443	ENDDO
	444	ELSEIF ( intermediate_timestep_count < &
	445	intermediate_timestep_count_max ) THEN
	446	DO k = nzb_v_inner(j,i)+1, nzt
[1337]	447	tv_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tv_m(k,j,i)
[736]	448	ENDDO
	449	ENDIF
	450	ENDIF
	451
	452	ENDIF
	453
	454	!
	455	!-- Tendency terms for w-velocity component
[1337]	456	tend(:,j,i) = 0.0_wp
[1001]	457	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	458	IF ( ws_scheme_mom ) THEN
	459	CALL advec_w_ws( i, j, i_omp_start, tn )
	460	ELSE
	461	CALL advec_w_pw( i, j )
	462	END IF
	463	ELSE
	464	CALL advec_w_up( i, j )
	465	ENDIF
[1001]	466	CALL diffusion_w( i, j )
[736]	467	CALL coriolis( i, j, 3 )
[940]	468
	469	IF ( .NOT. neutral ) THEN
	470	IF ( ocean ) THEN
[1179]	471	CALL buoyancy( i, j, rho, 3 )
[736]	472	ELSE
[940]	473	IF ( .NOT. humidity ) THEN
[1179]	474	CALL buoyancy( i, j, pt, 3 )
[940]	475	ELSE
[1179]	476	CALL buoyancy( i, j, vpt, 3 )
[940]	477	ENDIF
[736]	478	ENDIF
	479	ENDIF
	480
	481	!
	482	!-- Drag by plant canopy
	483	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 3 )
	484
	485	CALL user_actions( i, j, 'w-tendency' )
	486
	487	!
	488	!-- Prognostic equation for w-velocity component
	489	DO k = nzb_w_inner(j,i)+1, nzt-1
[1001]	490	w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	491	tsc(3) * tw_m(k,j,i) ) &
	492	- tsc(5) * rdf(k) * w(k,j,i)
[736]	493	ENDDO
	494
	495	!
	496	!-- Calculate tendencies for the next Runge-Kutta step
	497	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	498	IF ( intermediate_timestep_count == 1 ) THEN
	499	DO k = nzb_w_inner(j,i)+1, nzt-1
	500	tw_m(k,j,i) = tend(k,j,i)
	501	ENDDO
	502	ELSEIF ( intermediate_timestep_count < &
	503	intermediate_timestep_count_max ) THEN
	504	DO k = nzb_w_inner(j,i)+1, nzt-1
[1337]	505	tw_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tw_m(k,j,i)
[736]	506	ENDDO
	507	ENDIF
	508	ENDIF
[1361]	509
[736]	510	!
[940]	511	!-- If required, compute prognostic equation for potential temperature
	512	IF ( .NOT. neutral ) THEN
	513	!
	514	!-- Tendency terms for potential temperature
[1337]	515	tend(:,j,i) = 0.0_wp
[1001]	516	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[940]	517	IF ( ws_scheme_sca ) THEN
	518	CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, diss_s_pt, &
	519	flux_l_pt, diss_l_pt, i_omp_start, tn )
	520	ELSE
	521	CALL advec_s_pw( i, j, pt )
	522	ENDIF
	523	ELSE
	524	CALL advec_s_up( i, j, pt )
	525	ENDIF
[1001]	526	CALL diffusion_s( i, j, pt, shf, tswst, wall_heatflux )
[736]	527
	528	!
[940]	529	!-- If required compute heating/cooling due to long wave radiation
	530	!-- processes
	531	IF ( radiation ) THEN
	532	CALL calc_radiation( i, j )
	533	ENDIF
[736]	534
[1106]	535	!
[1361]	536	!-- If required compute impact of latent heat due to precipitation
	537	IF ( cloud_physics .AND. icloud_scheme == 1 .AND. &
	538	precipitation ) THEN
	539	CALL impact_of_latent_heat( i, j )
[940]	540	ENDIF
[736]	541
	542	!
[940]	543	!-- Consideration of heat sources within the plant canopy
[1337]	544	IF ( plant_canopy .AND. cthf /= 0.0_wp ) THEN
[940]	545	CALL plant_canopy_model( i, j, 4 )
	546	ENDIF
[736]	547
[940]	548	!
[1365]	549	!-- Large scale advection
	550	IF ( large_scale_forcing ) THEN
	551	CALL ls_advec( i, j, simulated_time, 'pt' )
	552	ENDIF
	553
	554	!
[1106]	555	!-- If required, compute effect of large-scale subsidence/ascent
[1365]	556	IF ( large_scale_subsidence .AND. &
	557	.NOT. use_subsidence_tendencies ) THEN
	558	CALL subsidence( i, j, tend, pt, pt_init, 2 )
[940]	559	ENDIF
[736]	560
[1241]	561	!
	562	!-- Nudging
	563	IF ( nudging ) CALL nudge( i, j, simulated_time, 'pt' )
	564
[940]	565	CALL user_actions( i, j, 'pt-tendency' )
[736]	566
	567	!
[940]	568	!-- Prognostic equation for potential temperature
	569	DO k = nzb_s_inner(j,i)+1, nzt
[1001]	570	pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	571	tsc(3) * tpt_m(k,j,i) ) &
	572	- tsc(5) * ( pt(k,j,i) - pt_init(k) ) *&
	573	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )
[940]	574	ENDDO
[736]	575
	576	!
[940]	577	!-- Calculate tendencies for the next Runge-Kutta step
	578	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	579	IF ( intermediate_timestep_count == 1 ) THEN
	580	DO k = nzb_s_inner(j,i)+1, nzt
	581	tpt_m(k,j,i) = tend(k,j,i)
	582	ENDDO
	583	ELSEIF ( intermediate_timestep_count < &
	584	intermediate_timestep_count_max ) THEN
	585	DO k = nzb_s_inner(j,i)+1, nzt
[1337]	586	tpt_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
	587	5.3125_wp * tpt_m(k,j,i)
[940]	588	ENDDO
	589	ENDIF
[736]	590	ENDIF
[940]	591
[736]	592	ENDIF
	593
	594	!
	595	!-- If required, compute prognostic equation for salinity
	596	IF ( ocean ) THEN
	597
	598	!
	599	!-- Tendency-terms for salinity
[1337]	600	tend(:,j,i) = 0.0_wp
[1001]	601	IF ( timestep_scheme(1:5) == 'runge' ) &
[736]	602	THEN
	603	IF ( ws_scheme_sca ) THEN
	604	CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa, &
	605	diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn )
	606	ELSE
	607	CALL advec_s_pw( i, j, sa )
	608	ENDIF
	609	ELSE
	610	CALL advec_s_up( i, j, sa )
	611	ENDIF
[1001]	612	CALL diffusion_s( i, j, sa, saswsb, saswst, wall_salinityflux )
[736]	613
	614	CALL user_actions( i, j, 'sa-tendency' )
	615
	616	!
	617	!-- Prognostic equation for salinity
	618	DO k = nzb_s_inner(j,i)+1, nzt
[1001]	619	sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	620	tsc(3) * tsa_m(k,j,i) ) &
	621	- tsc(5) * rdf_sc(k) * &
	622	( sa(k,j,i) - sa_init(k) )
[1337]	623	IF ( sa_p(k,j,i) < 0.0_wp ) sa_p(k,j,i) = 0.1_wp * sa(k,j,i)
[736]	624	ENDDO
	625
	626	!
	627	!-- Calculate tendencies for the next Runge-Kutta step
	628	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	629	IF ( intermediate_timestep_count == 1 ) THEN
	630	DO k = nzb_s_inner(j,i)+1, nzt
	631	tsa_m(k,j,i) = tend(k,j,i)
	632	ENDDO
	633	ELSEIF ( intermediate_timestep_count < &
	634	intermediate_timestep_count_max ) THEN
	635	DO k = nzb_s_inner(j,i)+1, nzt
[1337]	636	tsa_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
	637	5.3125_wp * tsa_m(k,j,i)
[736]	638	ENDDO
	639	ENDIF
	640	ENDIF
	641
	642	!
	643	!-- Calculate density by the equation of state for seawater
	644	CALL eqn_state_seawater( i, j )
	645
	646	ENDIF
	647
	648	!
	649	!-- If required, compute prognostic equation for total water content /
	650	!-- scalar
	651	IF ( humidity .OR. passive_scalar ) THEN
	652
	653	!
	654	!-- Tendency-terms for total water content / scalar
[1337]	655	tend(:,j,i) = 0.0_wp
[1001]	656	IF ( timestep_scheme(1:5) == 'runge' ) &
[736]	657	THEN
	658	IF ( ws_scheme_sca ) THEN
	659	CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
	660	diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
	661	ELSE
	662	CALL advec_s_pw( i, j, q )
	663	ENDIF
	664	ELSE
	665	CALL advec_s_up( i, j, q )
	666	ENDIF
[1001]	667	CALL diffusion_s( i, j, q, qsws, qswst, wall_qflux )
[1053]	668
[736]	669	!
[1361]	670	!-- If required compute decrease of total water content due to
	671	!-- precipitation
	672	IF ( cloud_physics .AND. icloud_scheme == 1 .AND. &
	673	precipitation ) THEN
	674	CALL calc_precipitation( i, j )
[736]	675	ENDIF
	676	!
	677	!-- Sink or source of scalar concentration due to canopy elements
[1106]	678	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 5 )
[736]	679
[1053]	680	!
[1365]	681	!-- Large scale advection
	682	IF ( large_scale_forcing ) THEN
	683	CALL ls_advec( i, j, simulated_time, 'q' )
	684	ENDIF
	685
	686	!
[736]	687	!-- If required compute influence of large-scale subsidence/ascent
[1365]	688	IF ( large_scale_subsidence .AND. &
	689	.NOT. use_subsidence_tendencies ) THEN
	690	CALL subsidence( i, j, tend, q, q_init, 3 )
[736]	691	ENDIF
	692
[1241]	693	!
	694	!-- Nudging
	695	IF ( nudging ) CALL nudge( i, j, simulated_time, 'q' )
	696
[736]	697	CALL user_actions( i, j, 'q-tendency' )
	698
	699	!
	700	!-- Prognostic equation for total water content / scalar
	701	DO k = nzb_s_inner(j,i)+1, nzt
[1001]	702	q_p(k,j,i) = q(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	703	tsc(3) * tq_m(k,j,i) ) &
	704	- tsc(5) * rdf_sc(k) * &
	705	( q(k,j,i) - q_init(k) )
[1337]	706	IF ( q_p(k,j,i) < 0.0_wp ) q_p(k,j,i) = 0.1_wp * q(k,j,i)
[736]	707	ENDDO
	708
	709	!
	710	!-- Calculate tendencies for the next Runge-Kutta step
	711	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	712	IF ( intermediate_timestep_count == 1 ) THEN
	713	DO k = nzb_s_inner(j,i)+1, nzt
	714	tq_m(k,j,i) = tend(k,j,i)
	715	ENDDO
	716	ELSEIF ( intermediate_timestep_count < &
	717	intermediate_timestep_count_max ) THEN
	718	DO k = nzb_s_inner(j,i)+1, nzt
[1337]	719	tq_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
	720	5.3125_wp * tq_m(k,j,i)
[736]	721	ENDDO
	722	ENDIF
	723	ENDIF
	724
[1053]	725	!
	726	!-- If required, calculate prognostic equations for rain water content
	727	!-- and rain drop concentration
[1115]	728	IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &
	729	precipitation ) THEN
[1053]	730	!
	731	!-- Calculate prognostic equation for rain water content
[1337]	732	tend(:,j,i) = 0.0_wp
[1053]	733	IF ( timestep_scheme(1:5) == 'runge' ) &
	734	THEN
	735	IF ( ws_scheme_sca ) THEN
	736	CALL advec_s_ws( i, j, qr, 'qr', flux_s_qr, &
	737	diss_s_qr, flux_l_qr, diss_l_qr, &
	738	i_omp_start, tn )
	739	ELSE
	740	CALL advec_s_pw( i, j, qr )
	741	ENDIF
	742	ELSE
	743	CALL advec_s_up( i, j, qr )
	744	ENDIF
	745	CALL diffusion_s( i, j, qr, qrsws, qrswst, wall_qrflux )
	746
[1115]	747	CALL user_actions( i, j, 'qr-tendency' )
[1053]	748	!
	749	!-- Prognostic equation for rain water content
	750	DO k = nzb_s_inner(j,i)+1, nzt
[1115]	751	qr_p(k,j,i) = qr(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	752	tsc(3) * tqr_m(k,j,i) ) &
	753	- tsc(5) * rdf_sc(k) * qr(k,j,i)
[1337]	754	IF ( qr_p(k,j,i) < 0.0_wp ) qr_p(k,j,i) = 0.0_wp
[1053]	755	ENDDO
	756	!
	757	!-- Calculate tendencies for the next Runge-Kutta step
	758	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	759	IF ( intermediate_timestep_count == 1 ) THEN
	760	DO k = nzb_s_inner(j,i)+1, nzt
	761	tqr_m(k,j,i) = tend(k,j,i)
	762	ENDDO
	763	ELSEIF ( intermediate_timestep_count < &
	764	intermediate_timestep_count_max ) THEN
	765	DO k = nzb_s_inner(j,i)+1, nzt
[1337]	766	tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
	767	5.3125_wp * tqr_m(k,j,i)
[1053]	768	ENDDO
	769	ENDIF
	770	ENDIF
	771
	772	!
	773	!-- Calculate prognostic equation for rain drop concentration.
[1337]	774	tend(:,j,i) = 0.0_wp
[1053]	775	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	776	IF ( ws_scheme_sca ) THEN
[1115]	777	CALL advec_s_ws( i, j, nr, 'nr', flux_s_nr, &
	778	diss_s_nr, flux_l_nr, diss_l_nr, &
	779	i_omp_start, tn )
[1053]	780	ELSE
	781	CALL advec_s_pw( i, j, nr )
	782	ENDIF
	783	ELSE
	784	CALL advec_s_up( i, j, nr )
	785	ENDIF
	786	CALL diffusion_s( i, j, nr, nrsws, nrswst, wall_nrflux )
	787
[1115]	788	CALL user_actions( i, j, 'nr-tendency' )
[1053]	789	!
	790	!-- Prognostic equation for rain drop concentration
	791	DO k = nzb_s_inner(j,i)+1, nzt
[1115]	792	nr_p(k,j,i) = nr(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	793	tsc(3) * tnr_m(k,j,i) ) &
	794	- tsc(5) * rdf_sc(k) * nr(k,j,i)
[1337]	795	IF ( nr_p(k,j,i) < 0.0_wp ) nr_p(k,j,i) = 0.0_wp
[1053]	796	ENDDO
	797	!
	798	!-- Calculate tendencies for the next Runge-Kutta step
	799	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	800	IF ( intermediate_timestep_count == 1 ) THEN
	801	DO k = nzb_s_inner(j,i)+1, nzt
	802	tnr_m(k,j,i) = tend(k,j,i)
	803	ENDDO
	804	ELSEIF ( intermediate_timestep_count < &
	805	intermediate_timestep_count_max ) THEN
	806	DO k = nzb_s_inner(j,i)+1, nzt
[1337]	807	tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
	808	5.3125_wp * tnr_m(k,j,i)
[1053]	809	ENDDO
	810	ENDIF
	811	ENDIF
	812
	813	ENDIF
	814
[1128]	815	ENDIF
	816
[736]	817	!
	818	!-- If required, compute prognostic equation for turbulent kinetic
	819	!-- energy (TKE)
	820	IF ( .NOT. constant_diffusion ) THEN
	821
	822	!
	823	!-- Tendency-terms for TKE
[1337]	824	tend(:,j,i) = 0.0_wp
[1332]	825	IF ( timestep_scheme(1:5) == 'runge' &
	826	.AND. .NOT. use_upstream_for_tke ) THEN
[736]	827	IF ( ws_scheme_sca ) THEN
[1001]	828	CALL advec_s_ws( i, j, e, 'e', flux_s_e, diss_s_e, &
	829	flux_l_e, diss_l_e , i_omp_start, tn )
[736]	830	ELSE
	831	CALL advec_s_pw( i, j, e )
	832	ENDIF
	833	ELSE
	834	CALL advec_s_up( i, j, e )
	835	ENDIF
[1001]	836	IF ( .NOT. humidity ) THEN
	837	IF ( ocean ) THEN
	838	CALL diffusion_e( i, j, prho, prho_reference )
[736]	839	ELSE
[1001]	840	CALL diffusion_e( i, j, pt, pt_reference )
[736]	841	ENDIF
	842	ELSE
[1001]	843	CALL diffusion_e( i, j, vpt, pt_reference )
[736]	844	ENDIF
	845	CALL production_e( i, j )
	846
	847	!
	848	!-- Additional sink term for flows through plant canopies
	849	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 6 )
	850
	851	CALL user_actions( i, j, 'e-tendency' )
	852
	853	!
	854	!-- Prognostic equation for TKE.
	855	!-- Eliminate negative TKE values, which can occur due to numerical
	856	!-- reasons in the course of the integration. In such cases the old
	857	!-- TKE value is reduced by 90%.
	858	DO k = nzb_s_inner(j,i)+1, nzt
[1001]	859	e_p(k,j,i) = e(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	860	tsc(3) * te_m(k,j,i) )
[1337]	861	IF ( e_p(k,j,i) < 0.0_wp ) e_p(k,j,i) = 0.1_wp * e(k,j,i)
[736]	862	ENDDO
	863
	864	!
	865	!-- Calculate tendencies for the next Runge-Kutta step
	866	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	867	IF ( intermediate_timestep_count == 1 ) THEN
	868	DO k = nzb_s_inner(j,i)+1, nzt
	869	te_m(k,j,i) = tend(k,j,i)
	870	ENDDO
	871	ELSEIF ( intermediate_timestep_count < &
	872	intermediate_timestep_count_max ) THEN
	873	DO k = nzb_s_inner(j,i)+1, nzt
[1337]	874	te_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
	875	5.3125_wp * te_m(k,j,i)
[736]	876	ENDDO
	877	ENDIF
	878	ENDIF
	879
	880	ENDIF ! TKE equation
	881
	882	ENDDO
	883	ENDDO
	884	!$OMP END PARALLEL
	885
	886	CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
	887
	888
	889	END SUBROUTINE prognostic_equations_cache
	890
	891
	892	SUBROUTINE prognostic_equations_vector
	893
	894	!------------------------------------------------------------------------------!
	895	! Version for vector machines
	896	!------------------------------------------------------------------------------!
	897
	898	IMPLICIT NONE
	899
[1320]	900	INTEGER(iwp) :: i !:
	901	INTEGER(iwp) :: j !:
	902	INTEGER(iwp) :: k !:
[736]	903
[1320]	904	REAL(wp) :: sbt !:
[736]	905
[1320]	906
[736]	907	!
[1361]	908	!-- If required, calculate cloud microphysical impacts (two-moment scheme)
	909	IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &
	910	( intermediate_timestep_count == 1 .OR. &
	911	call_microphysics_at_all_substeps ) &
	912	) THEN
	913	CALL cpu_log( log_point(51), 'microphysics', 'start' )
	914	CALL microphysics_control
	915	CALL cpu_log( log_point(51), 'microphysics', 'stop' )
	916	ENDIF
	917
	918	!
[736]	919	!-- u-velocity component
	920	CALL cpu_log( log_point(5), 'u-equation', 'start' )
	921
[1337]	922	tend = 0.0_wp
[1001]	923	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	924	IF ( ws_scheme_mom ) THEN
	925	CALL advec_u_ws
	926	ELSE
	927	CALL advec_u_pw
	928	ENDIF
	929	ELSE
[1001]	930	CALL advec_u_up
[736]	931	ENDIF
[1001]	932	CALL diffusion_u
[736]	933	CALL coriolis( 1 )
[940]	934	IF ( sloping_surface .AND. .NOT. neutral ) THEN
[1179]	935	CALL buoyancy( pt, 1 )
[940]	936	ENDIF
[736]	937
	938	!
	939	!-- Drag by plant canopy
	940	IF ( plant_canopy ) CALL plant_canopy_model( 1 )
	941
	942	!
	943	!-- External pressure gradient
	944	IF ( dp_external ) THEN
	945	DO i = nxlu, nxr
	946	DO j = nys, nyn
	947	DO k = dp_level_ind_b+1, nzt
	948	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
	949	ENDDO
	950	ENDDO
	951	ENDDO
	952	ENDIF
	953
[1241]	954	!
	955	!-- Nudging
	956	IF ( nudging ) CALL nudge( simulated_time, 'u' )
	957
[736]	958	CALL user_actions( 'u-tendency' )
	959
	960	!
	961	!-- Prognostic equation for u-velocity component
	962	DO i = nxlu, nxr
	963	DO j = nys, nyn
	964	DO k = nzb_u_inner(j,i)+1, nzt
[1001]	965	u_p(k,j,i) = u(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	966	tsc(3) * tu_m(k,j,i) ) &
	967	- tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
[736]	968	ENDDO
	969	ENDDO
	970	ENDDO
	971
	972	!
	973	!-- Calculate tendencies for the next Runge-Kutta step
	974	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	975	IF ( intermediate_timestep_count == 1 ) THEN
	976	DO i = nxlu, nxr
	977	DO j = nys, nyn
	978	DO k = nzb_u_inner(j,i)+1, nzt
	979	tu_m(k,j,i) = tend(k,j,i)
	980	ENDDO
	981	ENDDO
	982	ENDDO
	983	ELSEIF ( intermediate_timestep_count < &
	984	intermediate_timestep_count_max ) THEN
	985	DO i = nxlu, nxr
	986	DO j = nys, nyn
	987	DO k = nzb_u_inner(j,i)+1, nzt
[1337]	988	tu_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tu_m(k,j,i)
[736]	989	ENDDO
	990	ENDDO
	991	ENDDO
	992	ENDIF
	993	ENDIF
	994
	995	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
	996
	997	!
	998	!-- v-velocity component
	999	CALL cpu_log( log_point(6), 'v-equation', 'start' )
	1000
[1337]	1001	tend = 0.0_wp
[1001]	1002	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	1003	IF ( ws_scheme_mom ) THEN
	1004	CALL advec_v_ws
	1005	ELSE
	1006	CALL advec_v_pw
	1007	END IF
	1008	ELSE
[1001]	1009	CALL advec_v_up
[736]	1010	ENDIF
[1001]	1011	CALL diffusion_v
[736]	1012	CALL coriolis( 2 )
	1013
	1014	!
	1015	!-- Drag by plant canopy
	1016	IF ( plant_canopy ) CALL plant_canopy_model( 2 )
	1017
	1018	!
	1019	!-- External pressure gradient
	1020	IF ( dp_external ) THEN
	1021	DO i = nxl, nxr
	1022	DO j = nysv, nyn
	1023	DO k = dp_level_ind_b+1, nzt
	1024	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
	1025	ENDDO
	1026	ENDDO
	1027	ENDDO
	1028	ENDIF
	1029
[1241]	1030	!
	1031	!-- Nudging
	1032	IF ( nudging ) CALL nudge( simulated_time, 'v' )
	1033
[736]	1034	CALL user_actions( 'v-tendency' )
	1035
	1036	!
	1037	!-- Prognostic equation for v-velocity component
	1038	DO i = nxl, nxr
	1039	DO j = nysv, nyn
	1040	DO k = nzb_v_inner(j,i)+1, nzt
[1001]	1041	v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	1042	tsc(3) * tv_m(k,j,i) ) &
	1043	- tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
[736]	1044	ENDDO
	1045	ENDDO
	1046	ENDDO
	1047
	1048	!
	1049	!-- Calculate tendencies for the next Runge-Kutta step
	1050	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1051	IF ( intermediate_timestep_count == 1 ) THEN
	1052	DO i = nxl, nxr
	1053	DO j = nysv, nyn
	1054	DO k = nzb_v_inner(j,i)+1, nzt
	1055	tv_m(k,j,i) = tend(k,j,i)
	1056	ENDDO
	1057	ENDDO
	1058	ENDDO
	1059	ELSEIF ( intermediate_timestep_count < &
	1060	intermediate_timestep_count_max ) THEN
	1061	DO i = nxl, nxr
	1062	DO j = nysv, nyn
	1063	DO k = nzb_v_inner(j,i)+1, nzt
[1337]	1064	tv_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tv_m(k,j,i)
[736]	1065	ENDDO
	1066	ENDDO
	1067	ENDDO
	1068	ENDIF
	1069	ENDIF
	1070
	1071	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
	1072
	1073	!
	1074	!-- w-velocity component
	1075	CALL cpu_log( log_point(7), 'w-equation', 'start' )
	1076
[1353]	1077	tend = 0.0_wp
[1001]	1078	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	1079	IF ( ws_scheme_mom ) THEN
	1080	CALL advec_w_ws
	1081	ELSE
	1082	CALL advec_w_pw
	1083	ENDIF
	1084	ELSE
[1001]	1085	CALL advec_w_up
[736]	1086	ENDIF
[1001]	1087	CALL diffusion_w
[736]	1088	CALL coriolis( 3 )
[940]	1089
	1090	IF ( .NOT. neutral ) THEN
	1091	IF ( ocean ) THEN
[1179]	1092	CALL buoyancy( rho, 3 )
[736]	1093	ELSE
[940]	1094	IF ( .NOT. humidity ) THEN
[1179]	1095	CALL buoyancy( pt, 3 )
[940]	1096	ELSE
[1179]	1097	CALL buoyancy( vpt, 3 )
[940]	1098	ENDIF
[736]	1099	ENDIF
	1100	ENDIF
	1101
	1102	!
	1103	!-- Drag by plant canopy
	1104	IF ( plant_canopy ) CALL plant_canopy_model( 3 )
	1105
	1106	CALL user_actions( 'w-tendency' )
	1107
	1108	!
	1109	!-- Prognostic equation for w-velocity component
	1110	DO i = nxl, nxr
	1111	DO j = nys, nyn
	1112	DO k = nzb_w_inner(j,i)+1, nzt-1
[1001]	1113	w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	1114	tsc(3) * tw_m(k,j,i) ) &
	1115	- tsc(5) * rdf(k) * w(k,j,i)
[736]	1116	ENDDO
	1117	ENDDO
	1118	ENDDO
	1119
	1120	!
	1121	!-- Calculate tendencies for the next Runge-Kutta step
	1122	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1123	IF ( intermediate_timestep_count == 1 ) THEN
	1124	DO i = nxl, nxr
	1125	DO j = nys, nyn
	1126	DO k = nzb_w_inner(j,i)+1, nzt-1
	1127	tw_m(k,j,i) = tend(k,j,i)
	1128	ENDDO
	1129	ENDDO
	1130	ENDDO
	1131	ELSEIF ( intermediate_timestep_count < &
	1132	intermediate_timestep_count_max ) THEN
	1133	DO i = nxl, nxr
	1134	DO j = nys, nyn
	1135	DO k = nzb_w_inner(j,i)+1, nzt-1
[1337]	1136	tw_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tw_m(k,j,i)
[736]	1137	ENDDO
	1138	ENDDO
	1139	ENDDO
	1140	ENDIF
	1141	ENDIF
	1142
	1143	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
	1144
[940]	1145
[736]	1146	!
[940]	1147	!-- If required, compute prognostic equation for potential temperature
	1148	IF ( .NOT. neutral ) THEN
[736]	1149
[940]	1150	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
	1151
[736]	1152	!
[940]	1153	!-- pt-tendency terms with communication
	1154	sbt = tsc(2)
	1155	IF ( scalar_advec == 'bc-scheme' ) THEN
[736]	1156
[940]	1157	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
[736]	1158	!
[1001]	1159	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[1353]	1160	sbt = 1.0_wp
[940]	1161	ENDIF
[1337]	1162	tend = 0.0_wp
[940]	1163	CALL advec_s_bc( pt, 'pt' )
[1001]	1164
[736]	1165	ENDIF
[940]	1166
	1167	!
	1168	!-- pt-tendency terms with no communication
[1001]	1169	IF ( scalar_advec /= 'bc-scheme' ) THEN
[1337]	1170	tend = 0.0_wp
[1001]	1171	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[940]	1172	IF ( ws_scheme_sca ) THEN
	1173	CALL advec_s_ws( pt, 'pt' )
	1174	ELSE
	1175	CALL advec_s_pw( pt )
	1176	ENDIF
	1177	ELSE
[1001]	1178	CALL advec_s_up( pt )
[940]	1179	ENDIF
[736]	1180	ENDIF
	1181
[1001]	1182	CALL diffusion_s( pt, shf, tswst, wall_heatflux )
	1183
[736]	1184	!
[940]	1185	!-- If required compute heating/cooling due to long wave radiation processes
	1186	IF ( radiation ) THEN
	1187	CALL calc_radiation
	1188	ENDIF
[736]	1189
	1190	!
[940]	1191	!-- If required compute impact of latent heat due to precipitation
[1361]	1192	IF ( cloud_physics .AND. icloud_scheme == 1 .AND. precipitation ) THEN
[940]	1193	CALL impact_of_latent_heat
	1194	ENDIF
[736]	1195
	1196	!
[940]	1197	!-- Consideration of heat sources within the plant canopy
[1337]	1198	IF ( plant_canopy .AND. ( cthf /= 0.0_wp ) ) THEN
[940]	1199	CALL plant_canopy_model( 4 )
	1200	ENDIF
[736]	1201
[940]	1202	!
[1365]	1203	!-- Large scale advection
	1204	IF ( large_scale_forcing ) THEN
	1205	CALL ls_advec( simulated_time, 'pt' )
	1206	ENDIF
	1207
	1208	!
[940]	1209	!-- If required compute influence of large-scale subsidence/ascent
[1365]	1210	IF ( large_scale_subsidence .AND. &
	1211	.NOT. use_subsidence_tendencies ) THEN
	1212	CALL subsidence( tend, pt, pt_init, 2 )
[940]	1213	ENDIF
[736]	1214
[1241]	1215	!
	1216	!-- Nudging
	1217	IF ( nudging ) CALL nudge( simulated_time, 'pt' )
	1218
[940]	1219	CALL user_actions( 'pt-tendency' )
[736]	1220
	1221	!
[940]	1222	!-- Prognostic equation for potential temperature
	1223	DO i = nxl, nxr
	1224	DO j = nys, nyn
	1225	DO k = nzb_s_inner(j,i)+1, nzt
[1001]	1226	pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1227	tsc(3) * tpt_m(k,j,i) ) &
	1228	- tsc(5) * ( pt(k,j,i) - pt_init(k) ) *&
	1229	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )
[940]	1230	ENDDO
[736]	1231	ENDDO
	1232	ENDDO
	1233
	1234	!
[940]	1235	!-- Calculate tendencies for the next Runge-Kutta step
	1236	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1237	IF ( intermediate_timestep_count == 1 ) THEN
	1238	DO i = nxl, nxr
	1239	DO j = nys, nyn
	1240	DO k = nzb_s_inner(j,i)+1, nzt
	1241	tpt_m(k,j,i) = tend(k,j,i)
	1242	ENDDO
[736]	1243	ENDDO
	1244	ENDDO
[940]	1245	ELSEIF ( intermediate_timestep_count < &
	1246	intermediate_timestep_count_max ) THEN
	1247	DO i = nxl, nxr
	1248	DO j = nys, nyn
	1249	DO k = nzb_s_inner(j,i)+1, nzt
[1337]	1250	tpt_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
	1251	5.3125_wp * tpt_m(k,j,i)
[940]	1252	ENDDO
[736]	1253	ENDDO
	1254	ENDDO
[940]	1255	ENDIF
[736]	1256	ENDIF
[940]	1257
	1258	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
	1259
[736]	1260	ENDIF
	1261
	1262	!
	1263	!-- If required, compute prognostic equation for salinity
	1264	IF ( ocean ) THEN
	1265
	1266	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
	1267
	1268	!
	1269	!-- sa-tendency terms with communication
	1270	sbt = tsc(2)
	1271	IF ( scalar_advec == 'bc-scheme' ) THEN
	1272
	1273	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1274	!
[1001]	1275	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[1353]	1276	sbt = 1.0_wp
[736]	1277	ENDIF
[1337]	1278	tend = 0.0_wp
[736]	1279	CALL advec_s_bc( sa, 'sa' )
[1001]	1280
[736]	1281	ENDIF
	1282
	1283	!
	1284	!-- sa-tendency terms with no communication
[1001]	1285	IF ( scalar_advec /= 'bc-scheme' ) THEN
[1353]	1286	tend = 0.0_wp
[1001]	1287	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	1288	IF ( ws_scheme_sca ) THEN
	1289	CALL advec_s_ws( sa, 'sa' )
	1290	ELSE
	1291	CALL advec_s_pw( sa )
	1292	ENDIF
	1293	ELSE
[1001]	1294	CALL advec_s_up( sa )
[736]	1295	ENDIF
	1296	ENDIF
[1001]	1297
	1298	CALL diffusion_s( sa, saswsb, saswst, wall_salinityflux )
[736]	1299
	1300	CALL user_actions( 'sa-tendency' )
	1301
	1302	!
	1303	!-- Prognostic equation for salinity
	1304	DO i = nxl, nxr
	1305	DO j = nys, nyn
	1306	DO k = nzb_s_inner(j,i)+1, nzt
[1001]	1307	sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1308	tsc(3) * tsa_m(k,j,i) ) &
	1309	- tsc(5) * rdf_sc(k) * &
	1310	( sa(k,j,i) - sa_init(k) )
[1337]	1311	IF ( sa_p(k,j,i) < 0.0_wp ) sa_p(k,j,i) = 0.1_wp * sa(k,j,i)
[736]	1312	ENDDO
	1313	ENDDO
	1314	ENDDO
	1315
	1316	!
	1317	!-- Calculate tendencies for the next Runge-Kutta step
	1318	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1319	IF ( intermediate_timestep_count == 1 ) THEN
	1320	DO i = nxl, nxr
	1321	DO j = nys, nyn
	1322	DO k = nzb_s_inner(j,i)+1, nzt
	1323	tsa_m(k,j,i) = tend(k,j,i)
	1324	ENDDO
	1325	ENDDO
	1326	ENDDO
	1327	ELSEIF ( intermediate_timestep_count < &
	1328	intermediate_timestep_count_max ) THEN
	1329	DO i = nxl, nxr
	1330	DO j = nys, nyn
	1331	DO k = nzb_s_inner(j,i)+1, nzt
[1337]	1332	tsa_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
	1333	5.3125_wp * tsa_m(k,j,i)
[736]	1334	ENDDO
	1335	ENDDO
	1336	ENDDO
	1337	ENDIF
	1338	ENDIF
	1339
	1340	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
	1341
	1342	!
	1343	!-- Calculate density by the equation of state for seawater
	1344	CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
	1345	CALL eqn_state_seawater
	1346	CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
	1347
	1348	ENDIF
	1349
	1350	!
	1351	!-- If required, compute prognostic equation for total water content / scalar
	1352	IF ( humidity .OR. passive_scalar ) THEN
	1353
	1354	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
	1355
	1356	!
	1357	!-- Scalar/q-tendency terms with communication
	1358	sbt = tsc(2)
	1359	IF ( scalar_advec == 'bc-scheme' ) THEN
	1360
	1361	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1362	!
[1001]	1363	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[1337]	1364	sbt = 1.0_wp
[736]	1365	ENDIF
[1337]	1366	tend = 0.0_wp
[736]	1367	CALL advec_s_bc( q, 'q' )
[1001]	1368
[736]	1369	ENDIF
	1370
	1371	!
	1372	!-- Scalar/q-tendency terms with no communication
[1001]	1373	IF ( scalar_advec /= 'bc-scheme' ) THEN
[1337]	1374	tend = 0.0_wp
[1001]	1375	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	1376	IF ( ws_scheme_sca ) THEN
	1377	CALL advec_s_ws( q, 'q' )
	1378	ELSE
	1379	CALL advec_s_pw( q )
	1380	ENDIF
	1381	ELSE
[1001]	1382	CALL advec_s_up( q )
[736]	1383	ENDIF
	1384	ENDIF
[1001]	1385
	1386	CALL diffusion_s( q, qsws, qswst, wall_qflux )
[736]	1387
	1388	!
	1389	!-- If required compute decrease of total water content due to
	1390	!-- precipitation
[1361]	1391	IF ( cloud_physics .AND. icloud_scheme == 1 .AND. precipitation ) THEN
[736]	1392	CALL calc_precipitation
	1393	ENDIF
	1394
	1395	!
	1396	!-- Sink or source of scalar concentration due to canopy elements
	1397	IF ( plant_canopy ) CALL plant_canopy_model( 5 )
[1365]	1398
[736]	1399	!
[1365]	1400	!-- Large scale advection
	1401	IF ( large_scale_forcing ) THEN
	1402	CALL ls_advec( simulated_time, 'q' )
	1403	ENDIF
	1404
	1405	!
[736]	1406	!-- If required compute influence of large-scale subsidence/ascent
[1365]	1407	IF ( large_scale_subsidence .AND. &
	1408	.NOT. use_subsidence_tendencies ) THEN
	1409	CALL subsidence( tend, q, q_init, 3 )
[736]	1410	ENDIF
	1411
[1241]	1412	!
	1413	!-- Nudging
	1414	IF ( nudging ) CALL nudge( simulated_time, 'q' )
	1415
[736]	1416	CALL user_actions( 'q-tendency' )
	1417
	1418	!
	1419	!-- Prognostic equation for total water content / scalar
	1420	DO i = nxl, nxr
	1421	DO j = nys, nyn
	1422	DO k = nzb_s_inner(j,i)+1, nzt
[1001]	1423	q_p(k,j,i) = q(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1424	tsc(3) * tq_m(k,j,i) ) &
	1425	- tsc(5) * rdf_sc(k) * &
	1426	( q(k,j,i) - q_init(k) )
[1337]	1427	IF ( q_p(k,j,i) < 0.0_wp ) q_p(k,j,i) = 0.1_wp * q(k,j,i)
[736]	1428	ENDDO
	1429	ENDDO
	1430	ENDDO
	1431
	1432	!
	1433	!-- Calculate tendencies for the next Runge-Kutta step
	1434	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1435	IF ( intermediate_timestep_count == 1 ) THEN
	1436	DO i = nxl, nxr
	1437	DO j = nys, nyn
	1438	DO k = nzb_s_inner(j,i)+1, nzt
	1439	tq_m(k,j,i) = tend(k,j,i)
	1440	ENDDO
	1441	ENDDO
	1442	ENDDO
	1443	ELSEIF ( intermediate_timestep_count < &
	1444	intermediate_timestep_count_max ) THEN
	1445	DO i = nxl, nxr
	1446	DO j = nys, nyn
	1447	DO k = nzb_s_inner(j,i)+1, nzt
[1337]	1448	tq_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tq_m(k,j,i)
[736]	1449	ENDDO
	1450	ENDDO
	1451	ENDDO
	1452	ENDIF
	1453	ENDIF
	1454
	1455	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
	1456
[1361]	1457	!
	1458	!-- If required, calculate prognostic equations for rain water content
	1459	!-- and rain drop concentration
	1460	IF ( cloud_physics .AND. icloud_scheme == 0 .AND. precipitation ) THEN
	1461
	1462	CALL cpu_log( log_point(52), 'qr-equation', 'start' )
	1463
	1464	!
	1465	!-- Calculate prognostic equation for rain water content
	1466	sbt = tsc(2)
	1467	IF ( scalar_advec == 'bc-scheme' ) THEN
	1468
	1469	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1470	!
	1471	!-- Bott-Chlond scheme always uses Euler time step. Thus:
	1472	sbt = 1.0_wp
	1473	ENDIF
	1474	tend = 0.0_wp
	1475	CALL advec_s_bc( qr, 'qr' )
	1476
	1477	ENDIF
	1478
	1479	!
	1480	!-- qr-tendency terms with no communication
	1481	IF ( scalar_advec /= 'bc-scheme' ) THEN
	1482	tend = 0.0_wp
	1483	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1484	IF ( ws_scheme_sca ) THEN
	1485	CALL advec_s_ws( qr, 'qr' )
	1486	ELSE
	1487	CALL advec_s_pw( qr )
	1488	ENDIF
	1489	ELSE
	1490	CALL advec_s_up( qr )
	1491	ENDIF
	1492	ENDIF
	1493
	1494	CALL diffusion_s( qr, qrsws, qrswst, wall_qrflux )
	1495
	1496	CALL user_actions( 'qr-tendency' )
	1497
	1498	!
	1499	!-- Prognostic equation for rain water content
	1500	DO i = nxl, nxr
	1501	DO j = nys, nyn
	1502	DO k = nzb_s_inner(j,i)+1, nzt
	1503	qr_p(k,j,i) = qr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1504	tsc(3) * tqr_m(k,j,i) ) &
	1505	- tsc(5) * rdf_sc(k) * qr(k,j,i)
	1506	IF ( qr_p(k,j,i) < 0.0_wp ) qr_p(k,j,i) = 0.0_wp
	1507	ENDDO
	1508	ENDDO
	1509	ENDDO
	1510
	1511	!
	1512	!-- Calculate tendencies for the next Runge-Kutta step
	1513	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1514	IF ( intermediate_timestep_count == 1 ) THEN
	1515	DO i = nxl, nxr
	1516	DO j = nys, nyn
	1517	DO k = nzb_s_inner(j,i)+1, nzt
	1518	tqr_m(k,j,i) = tend(k,j,i)
	1519	ENDDO
	1520	ENDDO
	1521	ENDDO
	1522	ELSEIF ( intermediate_timestep_count < &
	1523	intermediate_timestep_count_max ) THEN
	1524	DO i = nxl, nxr
	1525	DO j = nys, nyn
	1526	DO k = nzb_s_inner(j,i)+1, nzt
	1527	tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * &
	1528	tqr_m(k,j,i)
	1529	ENDDO
	1530	ENDDO
	1531	ENDDO
	1532	ENDIF
	1533	ENDIF
	1534
	1535	CALL cpu_log( log_point(52), 'qr-equation', 'stop' )
	1536	CALL cpu_log( log_point(53), 'nr-equation', 'start' )
	1537
	1538	!
	1539	!-- Calculate prognostic equation for rain drop concentration
	1540	sbt = tsc(2)
	1541	IF ( scalar_advec == 'bc-scheme' ) THEN
	1542
	1543	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1544	!
	1545	!-- Bott-Chlond scheme always uses Euler time step. Thus:
	1546	sbt = 1.0_wp
	1547	ENDIF
	1548	tend = 0.0_wp
	1549	CALL advec_s_bc( nr, 'nr' )
	1550
	1551	ENDIF
	1552
	1553	!
	1554	!-- nr-tendency terms with no communication
	1555	IF ( scalar_advec /= 'bc-scheme' ) THEN
	1556	tend = 0.0_wp
	1557	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1558	IF ( ws_scheme_sca ) THEN
	1559	CALL advec_s_ws( nr, 'nr' )
	1560	ELSE
	1561	CALL advec_s_pw( nr )
	1562	ENDIF
	1563	ELSE
	1564	CALL advec_s_up( nr )
	1565	ENDIF
	1566	ENDIF
	1567
	1568	CALL diffusion_s( nr, nrsws, nrswst, wall_nrflux )
	1569
	1570	CALL user_actions( 'nr-tendency' )
	1571
	1572	!
	1573	!-- Prognostic equation for rain drop concentration
	1574	DO i = nxl, nxr
	1575	DO j = nys, nyn
	1576	DO k = nzb_s_inner(j,i)+1, nzt
	1577	nr_p(k,j,i) = nr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1578	tsc(3) * tnr_m(k,j,i) ) &
	1579	- tsc(5) * rdf_sc(k) * nr(k,j,i)
	1580	IF ( nr_p(k,j,i) < 0.0_wp ) nr_p(k,j,i) = 0.0_wp
	1581	ENDDO
	1582	ENDDO
	1583	ENDDO
	1584
	1585	!
	1586	!-- Calculate tendencies for the next Runge-Kutta step
	1587	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1588	IF ( intermediate_timestep_count == 1 ) THEN
	1589	DO i = nxl, nxr
	1590	DO j = nys, nyn
	1591	DO k = nzb_s_inner(j,i)+1, nzt
	1592	tnr_m(k,j,i) = tend(k,j,i)
	1593	ENDDO
	1594	ENDDO
	1595	ENDDO
	1596	ELSEIF ( intermediate_timestep_count < &
	1597	intermediate_timestep_count_max ) THEN
	1598	DO i = nxl, nxr
	1599	DO j = nys, nyn
	1600	DO k = nzb_s_inner(j,i)+1, nzt
	1601	tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * &
	1602	tnr_m(k,j,i)
	1603	ENDDO
	1604	ENDDO
	1605	ENDDO
	1606	ENDIF
	1607	ENDIF
	1608
	1609	CALL cpu_log( log_point(53), 'nr-equation', 'stop' )
	1610
	1611	ENDIF
	1612
[736]	1613	ENDIF
	1614
	1615	!
	1616	!-- If required, compute prognostic equation for turbulent kinetic
	1617	!-- energy (TKE)
	1618	IF ( .NOT. constant_diffusion ) THEN
	1619
	1620	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
	1621
	1622	sbt = tsc(2)
	1623	IF ( .NOT. use_upstream_for_tke ) THEN
	1624	IF ( scalar_advec == 'bc-scheme' ) THEN
	1625
	1626	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1627	!
[1001]	1628	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[1337]	1629	sbt = 1.0_wp
[736]	1630	ENDIF
[1337]	1631	tend = 0.0_wp
[736]	1632	CALL advec_s_bc( e, 'e' )
[1001]	1633
[736]	1634	ENDIF
	1635	ENDIF
	1636
	1637	!
	1638	!-- TKE-tendency terms with no communication
[1001]	1639	IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN
[736]	1640	IF ( use_upstream_for_tke ) THEN
[1337]	1641	tend = 0.0_wp
[736]	1642	CALL advec_s_up( e )
	1643	ELSE
[1337]	1644	tend = 0.0_wp
[1001]	1645	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	1646	IF ( ws_scheme_sca ) THEN
	1647	CALL advec_s_ws( e, 'e' )
	1648	ELSE
	1649	CALL advec_s_pw( e )
	1650	ENDIF
	1651	ELSE
[1001]	1652	CALL advec_s_up( e )
[736]	1653	ENDIF
	1654	ENDIF
[1001]	1655	ENDIF
	1656
	1657	IF ( .NOT. humidity ) THEN
	1658	IF ( ocean ) THEN
	1659	CALL diffusion_e( prho, prho_reference )
[736]	1660	ELSE
[1001]	1661	CALL diffusion_e( pt, pt_reference )
[736]	1662	ENDIF
[1001]	1663	ELSE
	1664	CALL diffusion_e( vpt, pt_reference )
[736]	1665	ENDIF
[1001]	1666
[736]	1667	CALL production_e
	1668
	1669	!
	1670	!-- Additional sink term for flows through plant canopies
	1671	IF ( plant_canopy ) CALL plant_canopy_model( 6 )
	1672	CALL user_actions( 'e-tendency' )
	1673
	1674	!
	1675	!-- Prognostic equation for TKE.
	1676	!-- Eliminate negative TKE values, which can occur due to numerical
	1677	!-- reasons in the course of the integration. In such cases the old TKE
	1678	!-- value is reduced by 90%.
	1679	DO i = nxl, nxr
	1680	DO j = nys, nyn
	1681	DO k = nzb_s_inner(j,i)+1, nzt
[1001]	1682	e_p(k,j,i) = e(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1683	tsc(3) * te_m(k,j,i) )
[1337]	1684	IF ( e_p(k,j,i) < 0.0_wp ) e_p(k,j,i) = 0.1_wp * e(k,j,i)
[736]	1685	ENDDO
	1686	ENDDO
	1687	ENDDO
	1688
	1689	!
	1690	!-- Calculate tendencies for the next Runge-Kutta step
	1691	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1692	IF ( intermediate_timestep_count == 1 ) THEN
	1693	DO i = nxl, nxr
	1694	DO j = nys, nyn
	1695	DO k = nzb_s_inner(j,i)+1, nzt
	1696	te_m(k,j,i) = tend(k,j,i)
	1697	ENDDO
	1698	ENDDO
	1699	ENDDO
	1700	ELSEIF ( intermediate_timestep_count < &
	1701	intermediate_timestep_count_max ) THEN
	1702	DO i = nxl, nxr
	1703	DO j = nys, nyn
	1704	DO k = nzb_s_inner(j,i)+1, nzt
[1337]	1705	te_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * te_m(k,j,i)
[736]	1706	ENDDO
	1707	ENDDO
	1708	ENDDO
	1709	ENDIF
	1710	ENDIF
	1711
	1712	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
	1713
	1714	ENDIF
	1715
	1716	END SUBROUTINE prognostic_equations_vector
	1717
	1718
[1015]	1719	SUBROUTINE prognostic_equations_acc
	1720
	1721	!------------------------------------------------------------------------------!
	1722	! Version for accelerator boards
	1723	!------------------------------------------------------------------------------!
	1724
	1725	IMPLICIT NONE
	1726
[1320]	1727	INTEGER(iwp) :: i !:
	1728	INTEGER(iwp) :: j !:
	1729	INTEGER(iwp) :: k !:
	1730	INTEGER(iwp) :: runge_step !:
[1015]	1731
[1320]	1732	REAL(wp) :: sbt !:
	1733
[1015]	1734	!
	1735	!-- Set switch for intermediate Runge-Kutta step
	1736	runge_step = 0
	1737	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1738	IF ( intermediate_timestep_count == 1 ) THEN
	1739	runge_step = 1
	1740	ELSEIF ( intermediate_timestep_count < &
	1741	intermediate_timestep_count_max ) THEN
	1742	runge_step = 2
	1743	ENDIF
	1744	ENDIF
	1745
	1746	!
[1361]	1747	!-- If required, calculate cloud microphysical impacts (two-moment scheme)
	1748	IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &
	1749	( intermediate_timestep_count == 1 .OR. &
	1750	call_microphysics_at_all_substeps ) &
	1751	) THEN
	1752	CALL cpu_log( log_point(51), 'microphysics', 'start' )
	1753	CALL microphysics_control
	1754	CALL cpu_log( log_point(51), 'microphysics', 'stop' )
	1755	ENDIF
	1756
	1757	!
[1015]	1758	!-- u-velocity component
	1759	!++ Statistics still not ported to accelerators
[1179]	1760	!$acc update device( hom, ref_state )
[1015]	1761	CALL cpu_log( log_point(5), 'u-equation', 'start' )
	1762
	1763	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1764	IF ( ws_scheme_mom ) THEN
	1765	CALL advec_u_ws_acc
	1766	ELSE
[1337]	1767	tend = 0.0_wp ! to be removed later??
[1015]	1768	CALL advec_u_pw
	1769	ENDIF
	1770	ELSE
	1771	CALL advec_u_up
	1772	ENDIF
	1773	CALL diffusion_u_acc
	1774	CALL coriolis_acc( 1 )
	1775	IF ( sloping_surface .AND. .NOT. neutral ) THEN
[1179]	1776	CALL buoyancy( pt, 1 )
[1015]	1777	ENDIF
	1778
	1779	!
	1780	!-- Drag by plant canopy
	1781	IF ( plant_canopy ) CALL plant_canopy_model( 1 )
	1782
	1783	!
	1784	!-- External pressure gradient
	1785	IF ( dp_external ) THEN
[1128]	1786	DO i = i_left, i_right
	1787	DO j = j_south, j_north
[1015]	1788	DO k = dp_level_ind_b+1, nzt
	1789	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
	1790	ENDDO
	1791	ENDDO
	1792	ENDDO
	1793	ENDIF
	1794
[1246]	1795	!
	1796	!-- Nudging
	1797	IF ( nudging ) CALL nudge( simulated_time, 'u' )
	1798
[1015]	1799	CALL user_actions( 'u-tendency' )
	1800
	1801	!
	1802	!-- Prognostic equation for u-velocity component
	1803	!$acc kernels present( nzb_u_inner, rdf, tend, tu_m, u, ug, u_p )
[1257]	1804	!$acc loop independent
[1128]	1805	DO i = i_left, i_right
[1257]	1806	!$acc loop independent
[1128]	1807	DO j = j_south, j_north
[1257]	1808	!$acc loop independent
[1015]	1809	DO k = 1, nzt
	1810	IF ( k > nzb_u_inner(j,i) ) THEN
	1811	u_p(k,j,i) = u(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	1812	tsc(3) * tu_m(k,j,i) ) &
	1813	- tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
	1814	!
	1815	!-- Tendencies for the next Runge-Kutta step
	1816	IF ( runge_step == 1 ) THEN
	1817	tu_m(k,j,i) = tend(k,j,i)
	1818	ELSEIF ( runge_step == 2 ) THEN
[1337]	1819	tu_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tu_m(k,j,i)
[1015]	1820	ENDIF
	1821	ENDIF
	1822	ENDDO
	1823	ENDDO
	1824	ENDDO
	1825	!$acc end kernels
	1826
	1827	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
	1828
	1829	!
	1830	!-- v-velocity component
	1831	CALL cpu_log( log_point(6), 'v-equation', 'start' )
	1832
	1833	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1834	IF ( ws_scheme_mom ) THEN
	1835	CALL advec_v_ws_acc
	1836	ELSE
[1337]	1837	tend = 0.0_wp ! to be removed later??
[1015]	1838	CALL advec_v_pw
	1839	END IF
	1840	ELSE
	1841	CALL advec_v_up
	1842	ENDIF
	1843	CALL diffusion_v_acc
	1844	CALL coriolis_acc( 2 )
	1845
	1846	!
	1847	!-- Drag by plant canopy
	1848	IF ( plant_canopy ) CALL plant_canopy_model( 2 )
	1849
	1850	!
	1851	!-- External pressure gradient
	1852	IF ( dp_external ) THEN
[1128]	1853	DO i = i_left, i_right
	1854	DO j = j_south, j_north
[1015]	1855	DO k = dp_level_ind_b+1, nzt
	1856	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
	1857	ENDDO
	1858	ENDDO
	1859	ENDDO
	1860	ENDIF
	1861
[1246]	1862	!
	1863	!-- Nudging
	1864	IF ( nudging ) CALL nudge( simulated_time, 'v' )
	1865
[1015]	1866	CALL user_actions( 'v-tendency' )
	1867
	1868	!
	1869	!-- Prognostic equation for v-velocity component
	1870	!$acc kernels present( nzb_v_inner, rdf, tend, tv_m, v, vg, v_p )
[1257]	1871	!$acc loop independent
[1128]	1872	DO i = i_left, i_right
[1257]	1873	!$acc loop independent
[1128]	1874	DO j = j_south, j_north
[1257]	1875	!$acc loop independent
[1015]	1876	DO k = 1, nzt
	1877	IF ( k > nzb_v_inner(j,i) ) THEN
	1878	v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	1879	tsc(3) * tv_m(k,j,i) ) &
	1880	- tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
	1881	!
	1882	!-- Tendencies for the next Runge-Kutta step
	1883	IF ( runge_step == 1 ) THEN
	1884	tv_m(k,j,i) = tend(k,j,i)
	1885	ELSEIF ( runge_step == 2 ) THEN
[1337]	1886	tv_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tv_m(k,j,i)
[1015]	1887	ENDIF
	1888	ENDIF
	1889	ENDDO
	1890	ENDDO
	1891	ENDDO
	1892	!$acc end kernels
	1893
	1894	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
	1895
	1896	!
	1897	!-- w-velocity component
	1898	CALL cpu_log( log_point(7), 'w-equation', 'start' )
	1899
	1900	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1901	IF ( ws_scheme_mom ) THEN
	1902	CALL advec_w_ws_acc
	1903	ELSE
[1337]	1904	tend = 0.0_wp ! to be removed later??
[1015]	1905	CALL advec_w_pw
	1906	ENDIF
	1907	ELSE
	1908	CALL advec_w_up
	1909	ENDIF
	1910	CALL diffusion_w_acc
	1911	CALL coriolis_acc( 3 )
	1912
	1913	IF ( .NOT. neutral ) THEN
	1914	IF ( ocean ) THEN
[1179]	1915	CALL buoyancy( rho, 3 )
[1015]	1916	ELSE
	1917	IF ( .NOT. humidity ) THEN
[1179]	1918	CALL buoyancy_acc( pt, 3 )
[1015]	1919	ELSE
[1179]	1920	CALL buoyancy( vpt, 3 )
[1015]	1921	ENDIF
	1922	ENDIF
	1923	ENDIF
	1924
	1925	!
	1926	!-- Drag by plant canopy
	1927	IF ( plant_canopy ) CALL plant_canopy_model( 3 )
	1928
	1929	CALL user_actions( 'w-tendency' )
	1930
	1931	!
	1932	!-- Prognostic equation for w-velocity component
	1933	!$acc kernels present( nzb_w_inner, rdf, tend, tw_m, w, w_p )
[1257]	1934	!$acc loop independent
[1128]	1935	DO i = i_left, i_right
[1257]	1936	!$acc loop independent
[1128]	1937	DO j = j_south, j_north
[1257]	1938	!$acc loop independent
[1015]	1939	DO k = 1, nzt-1
	1940	IF ( k > nzb_w_inner(j,i) ) THEN
	1941	w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	1942	tsc(3) * tw_m(k,j,i) ) &
	1943	- tsc(5) * rdf(k) * w(k,j,i)
	1944	!
	1945	!-- Tendencies for the next Runge-Kutta step
	1946	IF ( runge_step == 1 ) THEN
	1947	tw_m(k,j,i) = tend(k,j,i)
	1948	ELSEIF ( runge_step == 2 ) THEN
[1337]	1949	tw_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tw_m(k,j,i)
[1015]	1950	ENDIF
	1951	ENDIF
	1952	ENDDO
	1953	ENDDO
	1954	ENDDO
	1955	!$acc end kernels
	1956
	1957	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
	1958
	1959
	1960	!
	1961	!-- If required, compute prognostic equation for potential temperature
	1962	IF ( .NOT. neutral ) THEN
	1963
	1964	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
	1965
	1966	!
	1967	!-- pt-tendency terms with communication
	1968	sbt = tsc(2)
	1969	IF ( scalar_advec == 'bc-scheme' ) THEN
	1970
	1971	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1972	!
	1973	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[1353]	1974	sbt = 1.0_wp
[1015]	1975	ENDIF
[1337]	1976	tend = 0.0_wp
[1015]	1977	CALL advec_s_bc( pt, 'pt' )
	1978
	1979	ENDIF
	1980
	1981	!
	1982	!-- pt-tendency terms with no communication
	1983	IF ( scalar_advec /= 'bc-scheme' ) THEN
[1337]	1984	tend = 0.0_wp
[1015]	1985	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1986	IF ( ws_scheme_sca ) THEN
	1987	CALL advec_s_ws_acc( pt, 'pt' )
	1988	ELSE
[1337]	1989	tend = 0.0_wp ! to be removed later??
[1015]	1990	CALL advec_s_pw( pt )
	1991	ENDIF
	1992	ELSE
	1993	CALL advec_s_up( pt )
	1994	ENDIF
	1995	ENDIF
	1996
	1997	CALL diffusion_s_acc( pt, shf, tswst, wall_heatflux )
	1998
	1999	!
	2000	!-- If required compute heating/cooling due to long wave radiation processes
	2001	IF ( radiation ) THEN
	2002	CALL calc_radiation
	2003	ENDIF
	2004
	2005	!
	2006	!-- If required compute impact of latent heat due to precipitation
[1361]	2007	IF ( cloud_physics .AND. icloud_scheme == 1 .AND. precipitation ) THEN
[1015]	2008	CALL impact_of_latent_heat
	2009	ENDIF
	2010
	2011	!
	2012	!-- Consideration of heat sources within the plant canopy
[1337]	2013	IF ( plant_canopy .AND. ( cthf /= 0.0_wp ) ) THEN
[1015]	2014	CALL plant_canopy_model( 4 )
	2015	ENDIF
	2016
	2017	!
[1365]	2018	!-- Large scale advection
	2019	IF ( large_scale_forcing ) THEN
	2020	CALL ls_advec( simulated_time, 'pt' )
	2021	ENDIF
	2022
	2023	!
[1015]	2024	!-- If required compute influence of large-scale subsidence/ascent
[1365]	2025	IF ( large_scale_subsidence .AND. &
	2026	.NOT. use_subsidence_tendencies ) THEN
	2027	CALL subsidence( tend, pt, pt_init, 2 )
[1015]	2028	ENDIF
	2029
[1246]	2030	!
	2031	!-- Nudging
	2032	IF ( nudging ) CALL nudge( simulated_time, 'pt' )
	2033
[1015]	2034	CALL user_actions( 'pt-tendency' )
	2035
	2036	!
	2037	!-- Prognostic equation for potential temperature
	2038	!$acc kernels present( nzb_s_inner, rdf_sc, ptdf_x, ptdf_y, pt_init ) &
	2039	!$acc present( tend, tpt_m, pt, pt_p )
[1257]	2040	!$acc loop independent
[1128]	2041	DO i = i_left, i_right
[1257]	2042	!$acc loop independent
[1128]	2043	DO j = j_south, j_north
[1257]	2044	!$acc loop independent
[1015]	2045	DO k = 1, nzt
	2046	IF ( k > nzb_s_inner(j,i) ) THEN
	2047	pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	2048	tsc(3) * tpt_m(k,j,i) ) &
	2049	- tsc(5) * ( pt(k,j,i) - pt_init(k) ) *&
	2050	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )
	2051	!
	2052	!-- Tendencies for the next Runge-Kutta step
	2053	IF ( runge_step == 1 ) THEN
	2054	tpt_m(k,j,i) = tend(k,j,i)
	2055	ELSEIF ( runge_step == 2 ) THEN
[1337]	2056	tpt_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tpt_m(k,j,i)
[1015]	2057	ENDIF
	2058	ENDIF
	2059	ENDDO
	2060	ENDDO
	2061	ENDDO
	2062	!$acc end kernels
	2063
	2064	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
	2065
	2066	ENDIF
	2067
	2068	!
	2069	!-- If required, compute prognostic equation for salinity
	2070	IF ( ocean ) THEN
	2071
	2072	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
	2073
	2074	!
	2075	!-- sa-tendency terms with communication
	2076	sbt = tsc(2)
	2077	IF ( scalar_advec == 'bc-scheme' ) THEN
	2078
	2079	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	2080	!
	2081	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[1337]	2082	sbt = 1.0_wp
[1015]	2083	ENDIF
[1337]	2084	tend = 0.0_wp
[1015]	2085	CALL advec_s_bc( sa, 'sa' )
	2086
	2087	ENDIF
	2088
	2089	!
	2090	!-- sa-tendency terms with no communication
	2091	IF ( scalar_advec /= 'bc-scheme' ) THEN
[1337]	2092	tend = 0.0_wp
[1015]	2093	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	2094	IF ( ws_scheme_sca ) THEN
	2095	CALL advec_s_ws( sa, 'sa' )
	2096	ELSE
	2097	CALL advec_s_pw( sa )
	2098	ENDIF
	2099	ELSE
	2100	CALL advec_s_up( sa )
	2101	ENDIF
	2102	ENDIF
	2103
	2104	CALL diffusion_s( sa, saswsb, saswst, wall_salinityflux )
	2105
	2106	CALL user_actions( 'sa-tendency' )
	2107
	2108	!
	2109	!-- Prognostic equation for salinity
[1128]	2110	DO i = i_left, i_right
	2111	DO j = j_south, j_north
[1015]	2112	DO k = nzb_s_inner(j,i)+1, nzt
	2113	sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	2114	tsc(3) * tsa_m(k,j,i) ) &
	2115	- tsc(5) * rdf_sc(k) * &
	2116	( sa(k,j,i) - sa_init(k) )
[1337]	2117	IF ( sa_p(k,j,i) < 0.0_wp ) sa_p(k,j,i) = 0.1_wp * sa(k,j,i)
[1015]	2118	!
	2119	!-- Tendencies for the next Runge-Kutta step
	2120	IF ( runge_step == 1 ) THEN
	2121	tsa_m(k,j,i) = tend(k,j,i)
	2122	ELSEIF ( runge_step == 2 ) THEN
[1337]	2123	tsa_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tsa_m(k,j,i)
[1015]	2124	ENDIF
	2125	ENDDO
	2126	ENDDO
	2127	ENDDO
	2128
	2129	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
	2130
	2131	!
	2132	!-- Calculate density by the equation of state for seawater
	2133	CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
	2134	CALL eqn_state_seawater
	2135	CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
	2136
	2137	ENDIF
	2138
	2139	!
	2140	!-- If required, compute prognostic equation for total water content / scalar
	2141	IF ( humidity .OR. passive_scalar ) THEN
	2142
	2143	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
	2144
	2145	!
	2146	!-- Scalar/q-tendency terms with communication
	2147	sbt = tsc(2)
	2148	IF ( scalar_advec == 'bc-scheme' ) THEN
	2149
	2150	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	2151	!
	2152	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[1337]	2153	sbt = 1.0_wp
[1015]	2154	ENDIF
[1337]	2155	tend = 0.0_wp
[1015]	2156	CALL advec_s_bc( q, 'q' )
	2157
	2158	ENDIF
	2159
	2160	!
	2161	!-- Scalar/q-tendency terms with no communication
	2162	IF ( scalar_advec /= 'bc-scheme' ) THEN
[1337]	2163	tend = 0.0_wp
[1015]	2164	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	2165	IF ( ws_scheme_sca ) THEN
	2166	CALL advec_s_ws( q, 'q' )
	2167	ELSE
	2168	CALL advec_s_pw( q )
	2169	ENDIF
	2170	ELSE
	2171	CALL advec_s_up( q )
	2172	ENDIF
	2173	ENDIF
	2174
	2175	CALL diffusion_s( q, qsws, qswst, wall_qflux )
	2176
	2177	!
	2178	!-- If required compute decrease of total water content due to
	2179	!-- precipitation
[1361]	2180	IF ( cloud_physics .AND. icloud_scheme == 1 .AND. precipitation ) THEN
[1015]	2181	CALL calc_precipitation
	2182	ENDIF
	2183
	2184	!
	2185	!-- Sink or source of scalar concentration due to canopy elements
	2186	IF ( plant_canopy ) CALL plant_canopy_model( 5 )
	2187
	2188	!
[1365]	2189	!-- Large scale advection
	2190	IF ( large_scale_forcing ) THEN
	2191	CALL ls_advec( simulated_time, 'q' )
	2192	ENDIF
	2193
	2194	!
[1015]	2195	!-- If required compute influence of large-scale subsidence/ascent
[1365]	2196	IF ( large_scale_subsidence .AND. &
	2197	.NOT. use_subsidence_tendencies ) THEN
	2198	CALL subsidence( tend, q, q_init, 3 )
[1015]	2199	ENDIF
	2200
[1246]	2201	!
	2202	!-- Nudging
	2203	IF ( nudging ) CALL nudge( simulated_time, 'q' )
	2204
[1015]	2205	CALL user_actions( 'q-tendency' )
	2206
	2207	!
	2208	!-- Prognostic equation for total water content / scalar
[1128]	2209	DO i = i_left, i_right
	2210	DO j = j_south, j_north
[1015]	2211	DO k = nzb_s_inner(j,i)+1, nzt
	2212	q_p(k,j,i) = q(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	2213	tsc(3) * tq_m(k,j,i) ) &
	2214	- tsc(5) * rdf_sc(k) * &
	2215	( q(k,j,i) - q_init(k) )
[1337]	2216	IF ( q_p(k,j,i) < 0.0_wp ) q_p(k,j,i) = 0.1_wp * q(k,j,i)
[1015]	2217	!
	2218	!-- Tendencies for the next Runge-Kutta step
	2219	IF ( runge_step == 1 ) THEN
	2220	tq_m(k,j,i) = tend(k,j,i)
	2221	ELSEIF ( runge_step == 2 ) THEN
[1337]	2222	tq_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tq_m(k,j,i)
[1015]	2223	ENDIF
	2224	ENDDO
	2225	ENDDO
	2226	ENDDO
	2227
	2228	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
	2229
[1361]	2230	!
	2231	!-- If required, calculate prognostic equations for rain water content
	2232	!-- and rain drop concentration
	2233	IF ( cloud_physics .AND. icloud_scheme == 0 .AND. precipitation ) THEN
	2234
	2235	CALL cpu_log( log_point(52), 'qr-equation', 'start' )
	2236	!
	2237	!-- qr-tendency terms with communication
	2238	sbt = tsc(2)
	2239	IF ( scalar_advec == 'bc-scheme' ) THEN
	2240
	2241	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	2242	!
	2243	!-- Bott-Chlond scheme always uses Euler time step. Thus:
	2244	sbt = 1.0_wp
	2245	ENDIF
	2246	tend = 0.0_wp
	2247	CALL advec_s_bc( qr, 'qr' )
	2248
	2249	ENDIF
	2250
	2251	!
	2252	!-- qr-tendency terms with no communication
	2253	IF ( scalar_advec /= 'bc-scheme' ) THEN
	2254	tend = 0.0_wp
	2255	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	2256	IF ( ws_scheme_sca ) THEN
	2257	CALL advec_s_ws( qr, 'qr' )
	2258	ELSE
	2259	CALL advec_s_pw( qr )
	2260	ENDIF
	2261	ELSE
	2262	CALL advec_s_up( qr )
	2263	ENDIF
	2264	ENDIF
	2265
	2266	CALL diffusion_s( qr, qrsws, qrswst, wall_qrflux )
	2267
	2268	CALL user_actions( 'qr-tendency' )
	2269
	2270	!
	2271	!-- Prognostic equation for rain water content
	2272	DO i = i_left, i_right
	2273	DO j = j_south, j_north
	2274	DO k = nzb_s_inner(j,i)+1, nzt
	2275	qr_p(k,j,i) = qr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	2276	tsc(3) * tqr_m(k,j,i) ) &
	2277	- tsc(5) * rdf_sc(k) * qr(k,j,i)
	2278	IF ( qr_p(k,j,i) < 0.0_wp ) qr_p(k,j,i) = 0.0_wp
	2279	!
	2280	!-- Tendencies for the next Runge-Kutta step
	2281	IF ( runge_step == 1 ) THEN
	2282	tqr_m(k,j,i) = tend(k,j,i)
	2283	ELSEIF ( runge_step == 2 ) THEN
	2284	tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * &
	2285	tqr_m(k,j,i)
	2286	ENDIF
	2287	ENDDO
	2288	ENDDO
	2289	ENDDO
	2290
	2291	CALL cpu_log( log_point(52), 'qr-equation', 'stop' )
	2292	CALL cpu_log( log_point(53), 'nr-equation', 'start' )
	2293
	2294	!
	2295	!-- nr-tendency terms with communication
	2296	sbt = tsc(2)
	2297	IF ( scalar_advec == 'bc-scheme' ) THEN
	2298
	2299	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	2300	!
	2301	!-- Bott-Chlond scheme always uses Euler time step. Thus:
	2302	sbt = 1.0_wp
	2303	ENDIF
	2304	tend = 0.0_wp
	2305	CALL advec_s_bc( nr, 'nr' )
	2306
	2307	ENDIF
	2308
	2309	!
	2310	!-- nr-tendency terms with no communication
	2311	IF ( scalar_advec /= 'bc-scheme' ) THEN
	2312	tend = 0.0_wp
	2313	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	2314	IF ( ws_scheme_sca ) THEN
	2315	CALL advec_s_ws( nr, 'nr' )
	2316	ELSE
	2317	CALL advec_s_pw( nr )
	2318	ENDIF
	2319	ELSE
	2320	CALL advec_s_up( nr )
	2321	ENDIF
	2322	ENDIF
	2323
	2324	CALL diffusion_s( nr, nrsws, nrswst, wall_nrflux )
	2325
	2326	CALL user_actions( 'nr-tendency' )
	2327
	2328	!
	2329	!-- Prognostic equation for rain drop concentration
	2330	DO i = i_left, i_right
	2331	DO j = j_south, j_north
	2332	DO k = nzb_s_inner(j,i)+1, nzt
	2333	nr_p(k,j,i) = nr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	2334	tsc(3) * tnr_m(k,j,i) ) &
	2335	- tsc(5) * rdf_sc(k) * nr(k,j,i)
	2336	IF ( nr_p(k,j,i) < 0.0_wp ) nr_p(k,j,i) = 0.0_wp
	2337	!
	2338	!-- Tendencies for the next Runge-Kutta step
	2339	IF ( runge_step == 1 ) THEN
	2340	tnr_m(k,j,i) = tend(k,j,i)
	2341	ELSEIF ( runge_step == 2 ) THEN
	2342	tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * &
	2343	tnr_m(k,j,i)
	2344	ENDIF
	2345	ENDDO
	2346	ENDDO
	2347	ENDDO
	2348
	2349	CALL cpu_log( log_point(53), 'nr-equation', 'stop' )
	2350
	2351	ENDIF
	2352
[1015]	2353	ENDIF
	2354
	2355	!
	2356	!-- If required, compute prognostic equation for turbulent kinetic
	2357	!-- energy (TKE)
	2358	IF ( .NOT. constant_diffusion ) THEN
	2359
	2360	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
	2361
	2362	sbt = tsc(2)
	2363	IF ( .NOT. use_upstream_for_tke ) THEN
	2364	IF ( scalar_advec == 'bc-scheme' ) THEN
	2365
	2366	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	2367	!
	2368	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[1337]	2369	sbt = 1.0_wp
[1015]	2370	ENDIF
[1337]	2371	tend = 0.0_wp
[1015]	2372	CALL advec_s_bc( e, 'e' )
	2373
	2374	ENDIF
	2375	ENDIF
	2376
	2377	!
	2378	!-- TKE-tendency terms with no communication
	2379	IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN
	2380	IF ( use_upstream_for_tke ) THEN
[1337]	2381	tend = 0.0_wp
[1015]	2382	CALL advec_s_up( e )
	2383	ELSE
	2384	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	2385	IF ( ws_scheme_sca ) THEN
	2386	CALL advec_s_ws_acc( e, 'e' )
	2387	ELSE
[1337]	2388	tend = 0.0_wp ! to be removed later??
[1015]	2389	CALL advec_s_pw( e )
	2390	ENDIF
	2391	ELSE
[1337]	2392	tend = 0.0_wp ! to be removed later??
[1015]	2393	CALL advec_s_up( e )
	2394	ENDIF
	2395	ENDIF
	2396	ENDIF
	2397
	2398	IF ( .NOT. humidity ) THEN
	2399	IF ( ocean ) THEN
	2400	CALL diffusion_e( prho, prho_reference )
	2401	ELSE
	2402	CALL diffusion_e_acc( pt, pt_reference )
	2403	ENDIF
	2404	ELSE
	2405	CALL diffusion_e( vpt, pt_reference )
	2406	ENDIF
	2407
	2408	CALL production_e_acc
	2409
	2410	!
	2411	!-- Additional sink term for flows through plant canopies
	2412	IF ( plant_canopy ) CALL plant_canopy_model( 6 )
	2413	CALL user_actions( 'e-tendency' )
	2414
	2415	!
	2416	!-- Prognostic equation for TKE.
	2417	!-- Eliminate negative TKE values, which can occur due to numerical
	2418	!-- reasons in the course of the integration. In such cases the old TKE
	2419	!-- value is reduced by 90%.
	2420	!$acc kernels present( e, e_p, nzb_s_inner, tend, te_m )
[1257]	2421	!$acc loop independent
[1128]	2422	DO i = i_left, i_right
[1257]	2423	!$acc loop independent
[1128]	2424	DO j = j_south, j_north
[1257]	2425	!$acc loop independent
[1015]	2426	DO k = 1, nzt
	2427	IF ( k > nzb_s_inner(j,i) ) THEN
	2428	e_p(k,j,i) = e(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	2429	tsc(3) * te_m(k,j,i) )
[1337]	2430	IF ( e_p(k,j,i) < 0.0_wp ) e_p(k,j,i) = 0.1_wp * e(k,j,i)
[1015]	2431	!
	2432	!-- Tendencies for the next Runge-Kutta step
	2433	IF ( runge_step == 1 ) THEN
	2434	te_m(k,j,i) = tend(k,j,i)
	2435	ELSEIF ( runge_step == 2 ) THEN
[1337]	2436	te_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * te_m(k,j,i)
[1015]	2437	ENDIF
	2438	ENDIF
	2439	ENDDO
	2440	ENDDO
	2441	ENDDO
	2442	!$acc end kernels
	2443
	2444	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
	2445
	2446	ENDIF
	2447
	2448	END SUBROUTINE prognostic_equations_acc
	2449
	2450
[736]	2451	END MODULE prognostic_equations_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |