Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |