Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |