Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |