Home

Context Navigation

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

Last change on this file since 1409 was 1409, checked in by suehring, 10 years ago
Bugfix: set inflow boundary conditions also if no humidity or passive_scalar is used
Property svn:keywords set to `Id`
File size: 83.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |