Home

Context Navigation

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

Last change on this file since 1826 was 1826, checked in by maronga, 8 years ago
further modularization of radiation model and plant canopy model
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	! ------------------
[1826]	21	! Renamed canopy model calls.
[1822]	22	!
[1485]	23	! Former revisions:
	24	! -----------------
	25	! $Id: prognostic_equations.f90 1826 2016-04-07 12:01:39Z maronga $
	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, &
[1826]	286	ONLY: cthf, plant_canopy, pcm_tendency
[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
[1826]	405	IF ( plant_canopy ) CALL pcm_tendency( i, j, 1 )
[736]	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
[1826]	464	IF ( plant_canopy ) CALL pcm_tendency( i, j, 2 )
[736]	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
[1826]	533	IF ( plant_canopy ) CALL pcm_tendency( i, j, 3 )
[736]	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
[1826]	588	CALL pcm_tendency( i, j, 4 )
[940]	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
[1826]	722	IF ( plant_canopy ) CALL pcm_tendency( 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
[1826]	890	IF ( plant_canopy ) CALL pcm_tendency( i, j, 6 )
[736]	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
[1826]	984	IF ( plant_canopy ) CALL pcm_tendency( 1 )
[736]	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
[1826]	1060	IF ( plant_canopy ) CALL pcm_tendency( 2 )
[736]	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
[1826]	1148	IF ( plant_canopy ) CALL pcm_tendency( 3 )
[736]	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
[1826]	1237	CALL pcm_tendency( 4 )
[940]	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
[1826]	1435	IF ( plant_canopy ) CALL pcm_tendency( 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
[1826]	1534	CALL user_actions( 'qr-tendency' )
	1535
[1361]	1536	!
	1537	!-- Prognostic equation for rain water content
	1538	DO i = nxl, nxr
	1539	DO j = nys, nyn
	1540	DO k = nzb_s_inner(j,i)+1, nzt
	1541	qr_p(k,j,i) = qr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1542	tsc(3) * tqr_m(k,j,i) ) &
	1543	- tsc(5) * rdf_sc(k) * qr(k,j,i)
	1544	IF ( qr_p(k,j,i) < 0.0_wp ) qr_p(k,j,i) = 0.0_wp
	1545	ENDDO
	1546	ENDDO
	1547	ENDDO
	1548
	1549	!
	1550	!-- Calculate tendencies for the next Runge-Kutta step
	1551	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1552	IF ( intermediate_timestep_count == 1 ) THEN
	1553	DO i = nxl, nxr
	1554	DO j = nys, nyn
	1555	DO k = nzb_s_inner(j,i)+1, nzt
	1556	tqr_m(k,j,i) = tend(k,j,i)
	1557	ENDDO
	1558	ENDDO
	1559	ENDDO
	1560	ELSEIF ( intermediate_timestep_count < &
	1561	intermediate_timestep_count_max ) THEN
	1562	DO i = nxl, nxr
	1563	DO j = nys, nyn
	1564	DO k = nzb_s_inner(j,i)+1, nzt
	1565	tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * &
	1566	tqr_m(k,j,i)
	1567	ENDDO
	1568	ENDDO
	1569	ENDDO
	1570	ENDIF
	1571	ENDIF
	1572
	1573	CALL cpu_log( log_point(52), 'qr-equation', 'stop' )
	1574	CALL cpu_log( log_point(53), 'nr-equation', 'start' )
	1575
	1576	!
	1577	!-- Calculate prognostic equation for rain drop concentration
	1578	sbt = tsc(2)
	1579	IF ( scalar_advec == 'bc-scheme' ) THEN
	1580
	1581	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1582	!
	1583	!-- Bott-Chlond scheme always uses Euler time step. Thus:
	1584	sbt = 1.0_wp
	1585	ENDIF
	1586	tend = 0.0_wp
	1587	CALL advec_s_bc( nr, 'nr' )
	1588
	1589	ENDIF
	1590
	1591	!
	1592	!-- nr-tendency terms with no communication
	1593	IF ( scalar_advec /= 'bc-scheme' ) THEN
	1594	tend = 0.0_wp
	1595	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1596	IF ( ws_scheme_sca ) THEN
	1597	CALL advec_s_ws( nr, 'nr' )
	1598	ELSE
	1599	CALL advec_s_pw( nr )
	1600	ENDIF
	1601	ELSE
	1602	CALL advec_s_up( nr )
	1603	ENDIF
	1604	ENDIF
	1605
	1606	CALL diffusion_s( nr, nrsws, nrswst, wall_nrflux )
	1607
	1608	!
	1609	!-- Prognostic equation for rain drop concentration
	1610	DO i = nxl, nxr
	1611	DO j = nys, nyn
	1612	DO k = nzb_s_inner(j,i)+1, nzt
	1613	nr_p(k,j,i) = nr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1614	tsc(3) * tnr_m(k,j,i) ) &
	1615	- tsc(5) * rdf_sc(k) * nr(k,j,i)
	1616	IF ( nr_p(k,j,i) < 0.0_wp ) nr_p(k,j,i) = 0.0_wp
	1617	ENDDO
	1618	ENDDO
	1619	ENDDO
	1620
	1621	!
	1622	!-- Calculate tendencies for the next Runge-Kutta step
	1623	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1624	IF ( intermediate_timestep_count == 1 ) THEN
	1625	DO i = nxl, nxr
	1626	DO j = nys, nyn
	1627	DO k = nzb_s_inner(j,i)+1, nzt
	1628	tnr_m(k,j,i) = tend(k,j,i)
	1629	ENDDO
	1630	ENDDO
	1631	ENDDO
	1632	ELSEIF ( intermediate_timestep_count < &
	1633	intermediate_timestep_count_max ) THEN
	1634	DO i = nxl, nxr
	1635	DO j = nys, nyn
	1636	DO k = nzb_s_inner(j,i)+1, nzt
	1637	tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * &
	1638	tnr_m(k,j,i)
	1639	ENDDO
	1640	ENDDO
	1641	ENDDO
	1642	ENDIF
	1643	ENDIF
	1644
	1645	CALL cpu_log( log_point(53), 'nr-equation', 'stop' )
	1646
	1647	ENDIF
	1648
[736]	1649	ENDIF
	1650
	1651	!
	1652	!-- If required, compute prognostic equation for turbulent kinetic
	1653	!-- energy (TKE)
	1654	IF ( .NOT. constant_diffusion ) THEN
	1655
	1656	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
	1657
	1658	sbt = tsc(2)
	1659	IF ( .NOT. use_upstream_for_tke ) THEN
	1660	IF ( scalar_advec == 'bc-scheme' ) THEN
	1661
	1662	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1663	!
[1001]	1664	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[1337]	1665	sbt = 1.0_wp
[736]	1666	ENDIF
[1337]	1667	tend = 0.0_wp
[736]	1668	CALL advec_s_bc( e, 'e' )
[1001]	1669
[736]	1670	ENDIF
	1671	ENDIF
	1672
	1673	!
	1674	!-- TKE-tendency terms with no communication
[1001]	1675	IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN
[736]	1676	IF ( use_upstream_for_tke ) THEN
[1337]	1677	tend = 0.0_wp
[736]	1678	CALL advec_s_up( e )
	1679	ELSE
[1337]	1680	tend = 0.0_wp
[1001]	1681	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	1682	IF ( ws_scheme_sca ) THEN
	1683	CALL advec_s_ws( e, 'e' )
	1684	ELSE
	1685	CALL advec_s_pw( e )
	1686	ENDIF
	1687	ELSE
[1001]	1688	CALL advec_s_up( e )
[736]	1689	ENDIF
	1690	ENDIF
[1001]	1691	ENDIF
	1692
	1693	IF ( .NOT. humidity ) THEN
	1694	IF ( ocean ) THEN
	1695	CALL diffusion_e( prho, prho_reference )
[736]	1696	ELSE
[1001]	1697	CALL diffusion_e( pt, pt_reference )
[736]	1698	ENDIF
[1001]	1699	ELSE
	1700	CALL diffusion_e( vpt, pt_reference )
[736]	1701	ENDIF
[1001]	1702
[736]	1703	CALL production_e
	1704
	1705	!
	1706	!-- Additional sink term for flows through plant canopies
[1826]	1707	IF ( plant_canopy ) CALL pcm_tendency( 6 )
[736]	1708	CALL user_actions( 'e-tendency' )
	1709
	1710	!
	1711	!-- Prognostic equation for TKE.
	1712	!-- Eliminate negative TKE values, which can occur due to numerical
	1713	!-- reasons in the course of the integration. In such cases the old TKE
	1714	!-- value is reduced by 90%.
	1715	DO i = nxl, nxr
	1716	DO j = nys, nyn
	1717	DO k = nzb_s_inner(j,i)+1, nzt
[1001]	1718	e_p(k,j,i) = e(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1719	tsc(3) * te_m(k,j,i) )
[1337]	1720	IF ( e_p(k,j,i) < 0.0_wp ) e_p(k,j,i) = 0.1_wp * e(k,j,i)
[736]	1721	ENDDO
	1722	ENDDO
	1723	ENDDO
	1724
	1725	!
	1726	!-- Calculate tendencies for the next Runge-Kutta step
	1727	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1728	IF ( intermediate_timestep_count == 1 ) THEN
	1729	DO i = nxl, nxr
	1730	DO j = nys, nyn
	1731	DO k = nzb_s_inner(j,i)+1, nzt
	1732	te_m(k,j,i) = tend(k,j,i)
	1733	ENDDO
	1734	ENDDO
	1735	ENDDO
	1736	ELSEIF ( intermediate_timestep_count < &
	1737	intermediate_timestep_count_max ) THEN
	1738	DO i = nxl, nxr
	1739	DO j = nys, nyn
	1740	DO k = nzb_s_inner(j,i)+1, nzt
[1337]	1741	te_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * te_m(k,j,i)
[736]	1742	ENDDO
	1743	ENDDO
	1744	ENDDO
	1745	ENDIF
	1746	ENDIF
	1747
	1748	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
	1749
	1750	ENDIF
	1751
	1752	END SUBROUTINE prognostic_equations_vector
	1753
	1754
[1015]	1755	!------------------------------------------------------------------------------!
[1691]	1756	! Description:
	1757	! ------------
[1682]	1758	!> Version for accelerator boards
[1015]	1759	!------------------------------------------------------------------------------!
[1691]	1760
[1682]	1761	SUBROUTINE prognostic_equations_acc
[1015]	1762
[1682]	1763
[1015]	1764	IMPLICIT NONE
	1765
[1682]	1766	INTEGER(iwp) :: i !<
	1767	INTEGER(iwp) :: j !<
	1768	INTEGER(iwp) :: k !<
	1769	INTEGER(iwp) :: runge_step !<
[1015]	1770
[1682]	1771	REAL(wp) :: sbt !<
[1320]	1772
[1015]	1773	!
	1774	!-- Set switch for intermediate Runge-Kutta step
	1775	runge_step = 0
	1776	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1777	IF ( intermediate_timestep_count == 1 ) THEN
	1778	runge_step = 1
	1779	ELSEIF ( intermediate_timestep_count < &
	1780	intermediate_timestep_count_max ) THEN
	1781	runge_step = 2
	1782	ENDIF
	1783	ENDIF
	1784
	1785	!
[1361]	1786	!-- If required, calculate cloud microphysical impacts (two-moment scheme)
[1822]	1787	IF ( cloud_physics .AND. .NOT. microphysics_sat_adjust .AND. &
[1361]	1788	( intermediate_timestep_count == 1 .OR. &
	1789	call_microphysics_at_all_substeps ) &
	1790	) THEN
	1791	CALL cpu_log( log_point(51), 'microphysics', 'start' )
	1792	CALL microphysics_control
	1793	CALL cpu_log( log_point(51), 'microphysics', 'stop' )
	1794	ENDIF
	1795
	1796	!
[1015]	1797	!-- u-velocity component
[1374]	1798	!++ Statistics still not completely ported to accelerators
[1179]	1799	!$acc update device( hom, ref_state )
[1015]	1800	CALL cpu_log( log_point(5), 'u-equation', 'start' )
	1801
	1802	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1803	IF ( ws_scheme_mom ) THEN
	1804	CALL advec_u_ws_acc
	1805	ELSE
[1337]	1806	tend = 0.0_wp ! to be removed later??
[1015]	1807	CALL advec_u_pw
	1808	ENDIF
	1809	ELSE
	1810	CALL advec_u_up
	1811	ENDIF
	1812	CALL diffusion_u_acc
	1813	CALL coriolis_acc( 1 )
	1814	IF ( sloping_surface .AND. .NOT. neutral ) THEN
[1179]	1815	CALL buoyancy( pt, 1 )
[1015]	1816	ENDIF
	1817
	1818	!
	1819	!-- Drag by plant canopy
[1826]	1820	IF ( plant_canopy ) CALL pcm_tendency( 1 )
[1015]	1821
	1822	!
	1823	!-- External pressure gradient
	1824	IF ( dp_external ) THEN
[1128]	1825	DO i = i_left, i_right
	1826	DO j = j_south, j_north
[1015]	1827	DO k = dp_level_ind_b+1, nzt
	1828	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
	1829	ENDDO
	1830	ENDDO
	1831	ENDDO
	1832	ENDIF
	1833
[1246]	1834	!
	1835	!-- Nudging
[1691]	1836	IF ( nudging ) CALL nudge( simulated_time, 'u' )
[1246]	1837
[1015]	1838	CALL user_actions( 'u-tendency' )
	1839
	1840	!
	1841	!-- Prognostic equation for u-velocity component
[1398]	1842	!$acc kernels present( nzb_u_inner, rdf, tend, tu_m, u, u_init, u_p )
[1257]	1843	!$acc loop independent
[1128]	1844	DO i = i_left, i_right
[1257]	1845	!$acc loop independent
[1128]	1846	DO j = j_south, j_north
[1257]	1847	!$acc loop independent
[1015]	1848	DO k = 1, nzt
	1849	IF ( k > nzb_u_inner(j,i) ) THEN
	1850	u_p(k,j,i) = u(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	1851	tsc(3) * tu_m(k,j,i) ) &
[1398]	1852	- tsc(5) * rdf(k) * ( u(k,j,i) - u_init(k) )
[1015]	1853	!
	1854	!-- Tendencies for the next Runge-Kutta step
	1855	IF ( runge_step == 1 ) THEN
	1856	tu_m(k,j,i) = tend(k,j,i)
	1857	ELSEIF ( runge_step == 2 ) THEN
[1337]	1858	tu_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tu_m(k,j,i)
[1015]	1859	ENDIF
	1860	ENDIF
	1861	ENDDO
	1862	ENDDO
	1863	ENDDO
	1864	!$acc end kernels
	1865
	1866	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
	1867
	1868	!
	1869	!-- v-velocity component
	1870	CALL cpu_log( log_point(6), 'v-equation', 'start' )
	1871
	1872	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1873	IF ( ws_scheme_mom ) THEN
	1874	CALL advec_v_ws_acc
	1875	ELSE
[1337]	1876	tend = 0.0_wp ! to be removed later??
[1015]	1877	CALL advec_v_pw
	1878	END IF
	1879	ELSE
	1880	CALL advec_v_up
	1881	ENDIF
	1882	CALL diffusion_v_acc
	1883	CALL coriolis_acc( 2 )
	1884
	1885	!
	1886	!-- Drag by plant canopy
[1826]	1887	IF ( plant_canopy ) CALL pcm_tendency( 2 )
[1015]	1888
	1889	!
	1890	!-- External pressure gradient
	1891	IF ( dp_external ) THEN
[1128]	1892	DO i = i_left, i_right
	1893	DO j = j_south, j_north
[1015]	1894	DO k = dp_level_ind_b+1, nzt
	1895	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
	1896	ENDDO
	1897	ENDDO
	1898	ENDDO
	1899	ENDIF
	1900
[1246]	1901	!
	1902	!-- Nudging
[1691]	1903	IF ( nudging ) CALL nudge( simulated_time, 'v' )
[1246]	1904
[1015]	1905	CALL user_actions( 'v-tendency' )
	1906
	1907	!
	1908	!-- Prognostic equation for v-velocity component
[1398]	1909	!$acc kernels present( nzb_v_inner, rdf, tend, tv_m, v, v_init, v_p )
[1257]	1910	!$acc loop independent
[1128]	1911	DO i = i_left, i_right
[1257]	1912	!$acc loop independent
[1128]	1913	DO j = j_south, j_north
[1257]	1914	!$acc loop independent
[1015]	1915	DO k = 1, nzt
	1916	IF ( k > nzb_v_inner(j,i) ) THEN
	1917	v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	1918	tsc(3) * tv_m(k,j,i) ) &
[1398]	1919	- tsc(5) * rdf(k) * ( v(k,j,i) - v_init(k) )
[1015]	1920	!
	1921	!-- Tendencies for the next Runge-Kutta step
	1922	IF ( runge_step == 1 ) THEN
	1923	tv_m(k,j,i) = tend(k,j,i)
	1924	ELSEIF ( runge_step == 2 ) THEN
[1337]	1925	tv_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tv_m(k,j,i)
[1015]	1926	ENDIF
	1927	ENDIF
	1928	ENDDO
	1929	ENDDO
	1930	ENDDO
	1931	!$acc end kernels
	1932
	1933	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
	1934
	1935	!
	1936	!-- w-velocity component
	1937	CALL cpu_log( log_point(7), 'w-equation', 'start' )
	1938
	1939	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1940	IF ( ws_scheme_mom ) THEN
	1941	CALL advec_w_ws_acc
	1942	ELSE
[1337]	1943	tend = 0.0_wp ! to be removed later??
[1015]	1944	CALL advec_w_pw
	1945	ENDIF
	1946	ELSE
	1947	CALL advec_w_up
	1948	ENDIF
	1949	CALL diffusion_w_acc
	1950	CALL coriolis_acc( 3 )
	1951
	1952	IF ( .NOT. neutral ) THEN
	1953	IF ( ocean ) THEN
[1179]	1954	CALL buoyancy( rho, 3 )
[1015]	1955	ELSE
	1956	IF ( .NOT. humidity ) THEN
[1179]	1957	CALL buoyancy_acc( pt, 3 )
[1015]	1958	ELSE
[1179]	1959	CALL buoyancy( vpt, 3 )
[1015]	1960	ENDIF
	1961	ENDIF
	1962	ENDIF
	1963
	1964	!
	1965	!-- Drag by plant canopy
[1826]	1966	IF ( plant_canopy ) CALL pcm_tendency( 3 )
[1015]	1967
	1968	CALL user_actions( 'w-tendency' )
	1969
	1970	!
	1971	!-- Prognostic equation for w-velocity component
	1972	!$acc kernels present( nzb_w_inner, rdf, tend, tw_m, w, w_p )
[1257]	1973	!$acc loop independent
[1128]	1974	DO i = i_left, i_right
[1257]	1975	!$acc loop independent
[1128]	1976	DO j = j_south, j_north
[1257]	1977	!$acc loop independent
[1015]	1978	DO k = 1, nzt-1
	1979	IF ( k > nzb_w_inner(j,i) ) THEN
	1980	w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	1981	tsc(3) * tw_m(k,j,i) ) &
	1982	- tsc(5) * rdf(k) * w(k,j,i)
	1983	!
	1984	!-- Tendencies for the next Runge-Kutta step
	1985	IF ( runge_step == 1 ) THEN
	1986	tw_m(k,j,i) = tend(k,j,i)
	1987	ELSEIF ( runge_step == 2 ) THEN
[1337]	1988	tw_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tw_m(k,j,i)
[1015]	1989	ENDIF
	1990	ENDIF
	1991	ENDDO
	1992	ENDDO
	1993	ENDDO
	1994	!$acc end kernels
	1995
	1996	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
	1997
	1998
	1999	!
	2000	!-- If required, compute prognostic equation for potential temperature
	2001	IF ( .NOT. neutral ) THEN
	2002
	2003	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
	2004
	2005	!
	2006	!-- pt-tendency terms with communication
	2007	sbt = tsc(2)
	2008	IF ( scalar_advec == 'bc-scheme' ) THEN
	2009
	2010	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	2011	!
	2012	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[1353]	2013	sbt = 1.0_wp
[1015]	2014	ENDIF
[1337]	2015	tend = 0.0_wp
[1015]	2016	CALL advec_s_bc( pt, 'pt' )
	2017
	2018	ENDIF
	2019
	2020	!
	2021	!-- pt-tendency terms with no communication
	2022	IF ( scalar_advec /= 'bc-scheme' ) THEN
[1337]	2023	tend = 0.0_wp
[1015]	2024	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	2025	IF ( ws_scheme_sca ) THEN
	2026	CALL advec_s_ws_acc( pt, 'pt' )
	2027	ELSE
[1337]	2028	tend = 0.0_wp ! to be removed later??
[1015]	2029	CALL advec_s_pw( pt )
	2030	ENDIF
	2031	ELSE
	2032	CALL advec_s_up( pt )
	2033	ENDIF
	2034	ENDIF
	2035
	2036	CALL diffusion_s_acc( pt, shf, tswst, wall_heatflux )
	2037
	2038	!
	2039	!-- If required compute heating/cooling due to long wave radiation processes
[1496]	2040	IF ( cloud_top_radiation ) THEN
[1015]	2041	CALL calc_radiation
	2042	ENDIF
	2043
	2044	!
	2045	!-- Consideration of heat sources within the plant canopy
[1337]	2046	IF ( plant_canopy .AND. ( cthf /= 0.0_wp ) ) THEN
[1826]	2047	CALL pcm_tendency( 4 )
[1015]	2048	ENDIF
	2049
	2050	!
[1365]	2051	!-- Large scale advection
	2052	IF ( large_scale_forcing ) THEN
	2053	CALL ls_advec( simulated_time, 'pt' )
	2054	ENDIF
	2055
	2056	!
[1380]	2057	!-- Nudging
	2058	IF ( nudging ) CALL nudge( simulated_time, 'pt' )
	2059
	2060	!
[1015]	2061	!-- If required compute influence of large-scale subsidence/ascent
[1365]	2062	IF ( large_scale_subsidence .AND. &
	2063	.NOT. use_subsidence_tendencies ) THEN
	2064	CALL subsidence( tend, pt, pt_init, 2 )
[1015]	2065	ENDIF
	2066
[1691]	2067	IF ( radiation .AND. radiation_scheme == 'rrtmg' .AND. &
	2068	simulated_time > skip_time_do_radiation ) THEN
	2069	CALL radiation_tendency ( tend )
	2070	ENDIF
	2071
[1015]	2072	CALL user_actions( 'pt-tendency' )
	2073
	2074	!
	2075	!-- Prognostic equation for potential temperature
	2076	!$acc kernels present( nzb_s_inner, rdf_sc, ptdf_x, ptdf_y, pt_init ) &
	2077	!$acc present( tend, tpt_m, pt, pt_p )
[1257]	2078	!$acc loop independent
[1128]	2079	DO i = i_left, i_right
[1257]	2080	!$acc loop independent
[1128]	2081	DO j = j_south, j_north
[1257]	2082	!$acc loop independent
[1015]	2083	DO k = 1, nzt
	2084	IF ( k > nzb_s_inner(j,i) ) THEN
	2085	pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	2086	tsc(3) * tpt_m(k,j,i) ) &
	2087	- tsc(5) * ( pt(k,j,i) - pt_init(k) ) *&
	2088	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )
	2089	!
	2090	!-- Tendencies for the next Runge-Kutta step
	2091	IF ( runge_step == 1 ) THEN
	2092	tpt_m(k,j,i) = tend(k,j,i)
	2093	ELSEIF ( runge_step == 2 ) THEN
[1337]	2094	tpt_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tpt_m(k,j,i)
[1015]	2095	ENDIF
	2096	ENDIF
	2097	ENDDO
	2098	ENDDO
	2099	ENDDO
	2100	!$acc end kernels
	2101
	2102	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
	2103
	2104	ENDIF
	2105
	2106	!
	2107	!-- If required, compute prognostic equation for salinity
	2108	IF ( ocean ) THEN
	2109
	2110	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
	2111
	2112	!
	2113	!-- sa-tendency terms with communication
	2114	sbt = tsc(2)
	2115	IF ( scalar_advec == 'bc-scheme' ) THEN
	2116
	2117	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	2118	!
	2119	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[1337]	2120	sbt = 1.0_wp
[1015]	2121	ENDIF
[1337]	2122	tend = 0.0_wp
[1015]	2123	CALL advec_s_bc( sa, 'sa' )
	2124
	2125	ENDIF
	2126
	2127	!
	2128	!-- sa-tendency terms with no communication
	2129	IF ( scalar_advec /= 'bc-scheme' ) THEN
[1337]	2130	tend = 0.0_wp
[1015]	2131	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	2132	IF ( ws_scheme_sca ) THEN
	2133	CALL advec_s_ws( sa, 'sa' )
	2134	ELSE
	2135	CALL advec_s_pw( sa )
	2136	ENDIF
	2137	ELSE
	2138	CALL advec_s_up( sa )
	2139	ENDIF
	2140	ENDIF
	2141
	2142	CALL diffusion_s( sa, saswsb, saswst, wall_salinityflux )
	2143
	2144	CALL user_actions( 'sa-tendency' )
	2145
	2146	!
	2147	!-- Prognostic equation for salinity
[1128]	2148	DO i = i_left, i_right
	2149	DO j = j_south, j_north
[1015]	2150	DO k = nzb_s_inner(j,i)+1, nzt
	2151	sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	2152	tsc(3) * tsa_m(k,j,i) ) &
	2153	- tsc(5) * rdf_sc(k) * &
	2154	( sa(k,j,i) - sa_init(k) )
[1337]	2155	IF ( sa_p(k,j,i) < 0.0_wp ) sa_p(k,j,i) = 0.1_wp * sa(k,j,i)
[1015]	2156	!
	2157	!-- Tendencies for the next Runge-Kutta step
	2158	IF ( runge_step == 1 ) THEN
	2159	tsa_m(k,j,i) = tend(k,j,i)
	2160	ELSEIF ( runge_step == 2 ) THEN
[1337]	2161	tsa_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tsa_m(k,j,i)
[1015]	2162	ENDIF
	2163	ENDDO
	2164	ENDDO
	2165	ENDDO
	2166
	2167	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
	2168
	2169	!
	2170	!-- Calculate density by the equation of state for seawater
	2171	CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
	2172	CALL eqn_state_seawater
	2173	CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
	2174
	2175	ENDIF
	2176
	2177	!
	2178	!-- If required, compute prognostic equation for total water content / scalar
	2179	IF ( humidity .OR. passive_scalar ) THEN
	2180
	2181	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
	2182
	2183	!
	2184	!-- Scalar/q-tendency terms with communication
	2185	sbt = tsc(2)
	2186	IF ( scalar_advec == 'bc-scheme' ) THEN
	2187
	2188	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	2189	!
	2190	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[1337]	2191	sbt = 1.0_wp
[1015]	2192	ENDIF
[1337]	2193	tend = 0.0_wp
[1015]	2194	CALL advec_s_bc( q, 'q' )
	2195
	2196	ENDIF
	2197
	2198	!
	2199	!-- Scalar/q-tendency terms with no communication
	2200	IF ( scalar_advec /= 'bc-scheme' ) THEN
[1337]	2201	tend = 0.0_wp
[1015]	2202	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	2203	IF ( ws_scheme_sca ) THEN
	2204	CALL advec_s_ws( q, 'q' )
	2205	ELSE
	2206	CALL advec_s_pw( q )
	2207	ENDIF
	2208	ELSE
	2209	CALL advec_s_up( q )
	2210	ENDIF
	2211	ENDIF
	2212
	2213	CALL diffusion_s( q, qsws, qswst, wall_qflux )
	2214
	2215	!
	2216	!-- Sink or source of scalar concentration due to canopy elements
[1826]	2217	IF ( plant_canopy ) CALL pcm_tendency( 5 )
[1015]	2218
	2219	!
[1365]	2220	!-- Large scale advection
	2221	IF ( large_scale_forcing ) THEN
	2222	CALL ls_advec( simulated_time, 'q' )
	2223	ENDIF
	2224
	2225	!
[1380]	2226	!-- Nudging
	2227	IF ( nudging ) CALL nudge( simulated_time, 'q' )
	2228
	2229	!
[1015]	2230	!-- If required compute influence of large-scale subsidence/ascent
[1365]	2231	IF ( large_scale_subsidence .AND. &
	2232	.NOT. use_subsidence_tendencies ) THEN
	2233	CALL subsidence( tend, q, q_init, 3 )
[1015]	2234	ENDIF
	2235
	2236	CALL user_actions( 'q-tendency' )
	2237
	2238	!
	2239	!-- Prognostic equation for total water content / scalar
[1128]	2240	DO i = i_left, i_right
	2241	DO j = j_south, j_north
[1015]	2242	DO k = nzb_s_inner(j,i)+1, nzt
	2243	q_p(k,j,i) = q(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	2244	tsc(3) * tq_m(k,j,i) ) &
	2245	- tsc(5) * rdf_sc(k) * &
	2246	( q(k,j,i) - q_init(k) )
[1337]	2247	IF ( q_p(k,j,i) < 0.0_wp ) q_p(k,j,i) = 0.1_wp * q(k,j,i)
[1015]	2248	!
	2249	!-- Tendencies for the next Runge-Kutta step
	2250	IF ( runge_step == 1 ) THEN
	2251	tq_m(k,j,i) = tend(k,j,i)
	2252	ELSEIF ( runge_step == 2 ) THEN
[1337]	2253	tq_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tq_m(k,j,i)
[1015]	2254	ENDIF
	2255	ENDDO
	2256	ENDDO
	2257	ENDDO
	2258
	2259	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
	2260
[1361]	2261	!
	2262	!-- If required, calculate prognostic equations for rain water content
	2263	!-- and rain drop concentration
[1822]	2264	IF ( cloud_physics .AND. microphysics_seifert ) THEN
[1361]	2265
	2266	CALL cpu_log( log_point(52), 'qr-equation', 'start' )
	2267	!
	2268	!-- qr-tendency terms with communication
	2269	sbt = tsc(2)
	2270	IF ( scalar_advec == 'bc-scheme' ) THEN
	2271
	2272	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	2273	!
	2274	!-- Bott-Chlond scheme always uses Euler time step. Thus:
	2275	sbt = 1.0_wp
	2276	ENDIF
	2277	tend = 0.0_wp
	2278	CALL advec_s_bc( qr, 'qr' )
	2279
	2280	ENDIF
	2281
	2282	!
	2283	!-- qr-tendency terms with no communication
	2284	IF ( scalar_advec /= 'bc-scheme' ) THEN
	2285	tend = 0.0_wp
	2286	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	2287	IF ( ws_scheme_sca ) THEN
	2288	CALL advec_s_ws( qr, 'qr' )
	2289	ELSE
	2290	CALL advec_s_pw( qr )
	2291	ENDIF
	2292	ELSE
	2293	CALL advec_s_up( qr )
	2294	ENDIF
	2295	ENDIF
	2296
	2297	CALL diffusion_s( qr, qrsws, qrswst, wall_qrflux )
	2298
	2299	!
	2300	!-- Prognostic equation for rain water content
	2301	DO i = i_left, i_right
	2302	DO j = j_south, j_north
	2303	DO k = nzb_s_inner(j,i)+1, nzt
	2304	qr_p(k,j,i) = qr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	2305	tsc(3) * tqr_m(k,j,i) ) &
	2306	- tsc(5) * rdf_sc(k) * qr(k,j,i)
	2307	IF ( qr_p(k,j,i) < 0.0_wp ) qr_p(k,j,i) = 0.0_wp
	2308	!
	2309	!-- Tendencies for the next Runge-Kutta step
	2310	IF ( runge_step == 1 ) THEN
	2311	tqr_m(k,j,i) = tend(k,j,i)
	2312	ELSEIF ( runge_step == 2 ) THEN
	2313	tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * &
	2314	tqr_m(k,j,i)
	2315	ENDIF
	2316	ENDDO
	2317	ENDDO
	2318	ENDDO
	2319
	2320	CALL cpu_log( log_point(52), 'qr-equation', 'stop' )
	2321	CALL cpu_log( log_point(53), 'nr-equation', 'start' )
	2322
	2323	!
	2324	!-- nr-tendency terms with communication
	2325	sbt = tsc(2)
	2326	IF ( scalar_advec == 'bc-scheme' ) THEN
	2327
	2328	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	2329	!
	2330	!-- Bott-Chlond scheme always uses Euler time step. Thus:
	2331	sbt = 1.0_wp
	2332	ENDIF
	2333	tend = 0.0_wp
	2334	CALL advec_s_bc( nr, 'nr' )
	2335
	2336	ENDIF
	2337
	2338	!
	2339	!-- nr-tendency terms with no communication
	2340	IF ( scalar_advec /= 'bc-scheme' ) THEN
	2341	tend = 0.0_wp
	2342	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	2343	IF ( ws_scheme_sca ) THEN
	2344	CALL advec_s_ws( nr, 'nr' )
	2345	ELSE
	2346	CALL advec_s_pw( nr )
	2347	ENDIF
	2348	ELSE
	2349	CALL advec_s_up( nr )
	2350	ENDIF
	2351	ENDIF
	2352
	2353	CALL diffusion_s( nr, nrsws, nrswst, wall_nrflux )
	2354
	2355	!
	2356	!-- Prognostic equation for rain drop concentration
	2357	DO i = i_left, i_right
	2358	DO j = j_south, j_north
	2359	DO k = nzb_s_inner(j,i)+1, nzt
	2360	nr_p(k,j,i) = nr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	2361	tsc(3) * tnr_m(k,j,i) ) &
	2362	- tsc(5) * rdf_sc(k) * nr(k,j,i)
	2363	IF ( nr_p(k,j,i) < 0.0_wp ) nr_p(k,j,i) = 0.0_wp
	2364	!
	2365	!-- Tendencies for the next Runge-Kutta step
	2366	IF ( runge_step == 1 ) THEN
	2367	tnr_m(k,j,i) = tend(k,j,i)
	2368	ELSEIF ( runge_step == 2 ) THEN
	2369	tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * &
	2370	tnr_m(k,j,i)
	2371	ENDIF
	2372	ENDDO
	2373	ENDDO
	2374	ENDDO
	2375
	2376	CALL cpu_log( log_point(53), 'nr-equation', 'stop' )
	2377
	2378	ENDIF
	2379
[1015]	2380	ENDIF
	2381
	2382	!
	2383	!-- If required, compute prognostic equation for turbulent kinetic
	2384	!-- energy (TKE)
	2385	IF ( .NOT. constant_diffusion ) THEN
	2386
	2387	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
	2388
	2389	sbt = tsc(2)
	2390	IF ( .NOT. use_upstream_for_tke ) THEN
	2391	IF ( scalar_advec == 'bc-scheme' ) THEN
	2392
	2393	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	2394	!
	2395	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[1337]	2396	sbt = 1.0_wp
[1015]	2397	ENDIF
[1337]	2398	tend = 0.0_wp
[1015]	2399	CALL advec_s_bc( e, 'e' )
	2400
	2401	ENDIF
	2402	ENDIF
	2403
	2404	!
	2405	!-- TKE-tendency terms with no communication
	2406	IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN
	2407	IF ( use_upstream_for_tke ) THEN
[1337]	2408	tend = 0.0_wp
[1015]	2409	CALL advec_s_up( e )
	2410	ELSE
	2411	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	2412	IF ( ws_scheme_sca ) THEN
	2413	CALL advec_s_ws_acc( e, 'e' )
	2414	ELSE
[1337]	2415	tend = 0.0_wp ! to be removed later??
[1015]	2416	CALL advec_s_pw( e )
	2417	ENDIF
	2418	ELSE
[1337]	2419	tend = 0.0_wp ! to be removed later??
[1015]	2420	CALL advec_s_up( e )
	2421	ENDIF
	2422	ENDIF
	2423	ENDIF
	2424
	2425	IF ( .NOT. humidity ) THEN
	2426	IF ( ocean ) THEN
	2427	CALL diffusion_e( prho, prho_reference )
	2428	ELSE
	2429	CALL diffusion_e_acc( pt, pt_reference )
	2430	ENDIF
	2431	ELSE
	2432	CALL diffusion_e( vpt, pt_reference )
	2433	ENDIF
	2434
	2435	CALL production_e_acc
	2436
	2437	!
	2438	!-- Additional sink term for flows through plant canopies
[1826]	2439	IF ( plant_canopy ) CALL pcm_tendency( 6 )
[1015]	2440	CALL user_actions( 'e-tendency' )
	2441
	2442	!
	2443	!-- Prognostic equation for TKE.
	2444	!-- Eliminate negative TKE values, which can occur due to numerical
	2445	!-- reasons in the course of the integration. In such cases the old TKE
	2446	!-- value is reduced by 90%.
	2447	!$acc kernels present( e, e_p, nzb_s_inner, tend, te_m )
[1257]	2448	!$acc loop independent
[1128]	2449	DO i = i_left, i_right
[1257]	2450	!$acc loop independent
[1128]	2451	DO j = j_south, j_north
[1257]	2452	!$acc loop independent
[1015]	2453	DO k = 1, nzt
	2454	IF ( k > nzb_s_inner(j,i) ) THEN
	2455	e_p(k,j,i) = e(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	2456	tsc(3) * te_m(k,j,i) )
[1337]	2457	IF ( e_p(k,j,i) < 0.0_wp ) e_p(k,j,i) = 0.1_wp * e(k,j,i)
[1015]	2458	!
	2459	!-- Tendencies for the next Runge-Kutta step
	2460	IF ( runge_step == 1 ) THEN
	2461	te_m(k,j,i) = tend(k,j,i)
	2462	ELSEIF ( runge_step == 2 ) THEN
[1337]	2463	te_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * te_m(k,j,i)
[1015]	2464	ENDIF
	2465	ENDIF
	2466	ENDDO
	2467	ENDDO
	2468	ENDDO
	2469	!$acc end kernels
	2470
	2471	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
	2472
	2473	ENDIF
	2474
	2475	END SUBROUTINE prognostic_equations_acc
	2476
	2477
[736]	2478	END MODULE prognostic_equations_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |