Home

Context Navigation

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

Last change on this file since 2007 was 2007, checked in by kanani, 8 years ago
changes in the course of urban surface model implementation
Property svn:keywords set to `Id`
File size: 91.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |