Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |