Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |