Home

Context Navigation

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

Last change on this file since 1684 was 1683, checked in by knoop, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 84.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |