Home

Context Navigation

prognostic_equations.f90 @ 1484

Last change on this file since 1484 was 1484, checked in by kanani, 10 years ago

New:
---
Subroutine init_plant_canopy added to module plant_canopy_model_mod. (plant_canopy_model)
Alternative method for lad-profile construction added, also, new parameters added.
(header, package_parin, plant_canopy_model, read_var_list, write_var_list)
plant_canopy_model-dependency added to several subroutines. (Makefile)
New package/namelist canopy_par for canopy-related parameters added. (package_parin)

Changed:
---
Code structure of the plant canopy model changed, all canopy-model related code
combined to module plant_canopy_model_mod. (check_parameters, init_3d_model,
modules, timestep)
Module plant_canopy_model_mod added in USE-lists of some subroutines. (check_parameters,
header, init_3d_model, package_parin, read_var_list, user_init_plant_canopy, write_var_list)
Canopy initialization moved to new subroutine init_plant_canopy. (check_parameters,
init_3d_model, plant_canopy_model)
Calculation of canopy timestep-criterion removed, instead, the canopy
drag is now directly limited in the calculation of the canopy tendency terms.
(plant_canopy_model, timestep)
Some parameters renamed. (check_parameters, header, init_plant_canopy,
plant_canopy_model, read_var_list, write_var_list)
Unnecessary 3d-arrays removed. (init_plant_canopy, plant_canopy_model, user_init_plant_canopy)
Parameter checks regarding canopy initialization added. (check_parameters)
All canopy steering parameters moved from namelist inipar to canopy_par. (package_parin, parin)
Some redundant MPI communication removed. (init_plant_canopy)

Bugfix:
---
Missing KIND-attribute for REAL constant added. (check_parameters)
DO-WHILE-loop for lad-profile output restricted. (header)
Removed double-listing of use_upstream_for_tke in ONLY-list of module
control_parameters. (prognostic_equations)

Property svn:keywords set to Id

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

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

Context Navigation

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

Download in other formats: