Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |