Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |