Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |