Home

Context Navigation

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

Last change on this file since 1873 was 1873, checked in by maronga, 8 years ago
revised renaming of modules
Property svn:keywords set to `Id`
File size: 83.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |