Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |