Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |