Home

Context Navigation

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

Last change on this file since 1665 was 1586, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 84.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |