Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |