Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |