Home

Context Navigation

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

Last change on this file since 1332 was 1332, checked in by suehring, 11 years ago
Bugfix use_upstream_for_tke
Property svn:keywords set to `Id`
File size: 69.5 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	! ------------------
[1332]	22	! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke
	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: prognostic_equations.f90 1332 2014-03-25 11:59:43Z suehring $
	27	!
[1332]	28	! 1330 2014-03-24 17:29:32Z suehring
[1331]	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
[1332]	793	IF ( timestep_scheme(1:5) == 'runge' &
	794	.AND. .NOT. use_upstream_for_tke ) THEN
[736]	795	IF ( ws_scheme_sca ) THEN
[1001]	796	CALL advec_s_ws( i, j, e, 'e', flux_s_e, diss_s_e, &
	797	flux_l_e, diss_l_e , i_omp_start, tn )
[736]	798	ELSE
	799	CALL advec_s_pw( i, j, e )
	800	ENDIF
	801	ELSE
	802	CALL advec_s_up( i, j, e )
	803	ENDIF
[1001]	804	IF ( .NOT. humidity ) THEN
	805	IF ( ocean ) THEN
	806	CALL diffusion_e( i, j, prho, prho_reference )
[736]	807	ELSE
[1001]	808	CALL diffusion_e( i, j, pt, pt_reference )
[736]	809	ENDIF
	810	ELSE
[1001]	811	CALL diffusion_e( i, j, vpt, pt_reference )
[736]	812	ENDIF
	813	CALL production_e( i, j )
	814
	815	!
	816	!-- Additional sink term for flows through plant canopies
	817	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 6 )
	818
	819	CALL user_actions( i, j, 'e-tendency' )
	820
	821	!
	822	!-- Prognostic equation for TKE.
	823	!-- Eliminate negative TKE values, which can occur due to numerical
	824	!-- reasons in the course of the integration. In such cases the old
	825	!-- TKE value is reduced by 90%.
	826	DO k = nzb_s_inner(j,i)+1, nzt
[1001]	827	e_p(k,j,i) = e(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	828	tsc(3) * te_m(k,j,i) )
[736]	829	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
	830	ENDDO
	831
	832	!
	833	!-- Calculate tendencies for the next Runge-Kutta step
	834	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	835	IF ( intermediate_timestep_count == 1 ) THEN
	836	DO k = nzb_s_inner(j,i)+1, nzt
	837	te_m(k,j,i) = tend(k,j,i)
	838	ENDDO
	839	ELSEIF ( intermediate_timestep_count < &
	840	intermediate_timestep_count_max ) THEN
	841	DO k = nzb_s_inner(j,i)+1, nzt
	842	te_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	843	5.3125 * te_m(k,j,i)
	844	ENDDO
	845	ENDIF
	846	ENDIF
	847
	848	ENDIF ! TKE equation
	849
	850	ENDDO
	851	ENDDO
	852	!$OMP END PARALLEL
	853
	854	CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
	855
	856
	857	END SUBROUTINE prognostic_equations_cache
	858
	859
	860	SUBROUTINE prognostic_equations_vector
	861
	862	!------------------------------------------------------------------------------!
	863	! Version for vector machines
	864	!------------------------------------------------------------------------------!
	865
	866	IMPLICIT NONE
	867
[1320]	868	INTEGER(iwp) :: i !:
	869	INTEGER(iwp) :: j !:
	870	INTEGER(iwp) :: k !:
[736]	871
[1320]	872	REAL(wp) :: sbt !:
[736]	873
[1320]	874
[736]	875	!
	876	!-- u-velocity component
	877	CALL cpu_log( log_point(5), 'u-equation', 'start' )
	878
[1001]	879	tend = 0.0
	880	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	881	IF ( ws_scheme_mom ) THEN
	882	CALL advec_u_ws
	883	ELSE
	884	CALL advec_u_pw
	885	ENDIF
	886	ELSE
[1001]	887	CALL advec_u_up
[736]	888	ENDIF
[1001]	889	CALL diffusion_u
[736]	890	CALL coriolis( 1 )
[940]	891	IF ( sloping_surface .AND. .NOT. neutral ) THEN
[1179]	892	CALL buoyancy( pt, 1 )
[940]	893	ENDIF
[736]	894
	895	!
	896	!-- Drag by plant canopy
	897	IF ( plant_canopy ) CALL plant_canopy_model( 1 )
	898
	899	!
	900	!-- External pressure gradient
	901	IF ( dp_external ) THEN
	902	DO i = nxlu, nxr
	903	DO j = nys, nyn
	904	DO k = dp_level_ind_b+1, nzt
	905	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
	906	ENDDO
	907	ENDDO
	908	ENDDO
	909	ENDIF
	910
[1241]	911	!
	912	!-- Nudging
	913	IF ( nudging ) CALL nudge( simulated_time, 'u' )
	914
[736]	915	CALL user_actions( 'u-tendency' )
	916
	917	!
	918	!-- Prognostic equation for u-velocity component
	919	DO i = nxlu, nxr
	920	DO j = nys, nyn
	921	DO k = nzb_u_inner(j,i)+1, nzt
[1001]	922	u_p(k,j,i) = u(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	923	tsc(3) * tu_m(k,j,i) ) &
	924	- tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
[736]	925	ENDDO
	926	ENDDO
	927	ENDDO
	928
	929	!
	930	!-- Calculate tendencies for the next Runge-Kutta step
	931	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	932	IF ( intermediate_timestep_count == 1 ) THEN
	933	DO i = nxlu, nxr
	934	DO j = nys, nyn
	935	DO k = nzb_u_inner(j,i)+1, nzt
	936	tu_m(k,j,i) = tend(k,j,i)
	937	ENDDO
	938	ENDDO
	939	ENDDO
	940	ELSEIF ( intermediate_timestep_count < &
	941	intermediate_timestep_count_max ) THEN
	942	DO i = nxlu, nxr
	943	DO j = nys, nyn
	944	DO k = nzb_u_inner(j,i)+1, nzt
	945	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
	946	ENDDO
	947	ENDDO
	948	ENDDO
	949	ENDIF
	950	ENDIF
	951
	952	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
	953
	954	!
	955	!-- v-velocity component
	956	CALL cpu_log( log_point(6), 'v-equation', 'start' )
	957
[1001]	958	tend = 0.0
	959	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	960	IF ( ws_scheme_mom ) THEN
	961	CALL advec_v_ws
	962	ELSE
	963	CALL advec_v_pw
	964	END IF
	965	ELSE
[1001]	966	CALL advec_v_up
[736]	967	ENDIF
[1001]	968	CALL diffusion_v
[736]	969	CALL coriolis( 2 )
	970
	971	!
	972	!-- Drag by plant canopy
	973	IF ( plant_canopy ) CALL plant_canopy_model( 2 )
	974
	975	!
	976	!-- External pressure gradient
	977	IF ( dp_external ) THEN
	978	DO i = nxl, nxr
	979	DO j = nysv, nyn
	980	DO k = dp_level_ind_b+1, nzt
	981	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
	982	ENDDO
	983	ENDDO
	984	ENDDO
	985	ENDIF
	986
[1241]	987	!
	988	!-- Nudging
	989	IF ( nudging ) CALL nudge( simulated_time, 'v' )
	990
[736]	991	CALL user_actions( 'v-tendency' )
	992
	993	!
	994	!-- Prognostic equation for v-velocity component
	995	DO i = nxl, nxr
	996	DO j = nysv, nyn
	997	DO k = nzb_v_inner(j,i)+1, nzt
[1001]	998	v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	999	tsc(3) * tv_m(k,j,i) ) &
	1000	- tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
[736]	1001	ENDDO
	1002	ENDDO
	1003	ENDDO
	1004
	1005	!
	1006	!-- Calculate tendencies for the next Runge-Kutta step
	1007	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1008	IF ( intermediate_timestep_count == 1 ) THEN
	1009	DO i = nxl, nxr
	1010	DO j = nysv, nyn
	1011	DO k = nzb_v_inner(j,i)+1, nzt
	1012	tv_m(k,j,i) = tend(k,j,i)
	1013	ENDDO
	1014	ENDDO
	1015	ENDDO
	1016	ELSEIF ( intermediate_timestep_count < &
	1017	intermediate_timestep_count_max ) THEN
	1018	DO i = nxl, nxr
	1019	DO j = nysv, nyn
	1020	DO k = nzb_v_inner(j,i)+1, nzt
	1021	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
	1022	ENDDO
	1023	ENDDO
	1024	ENDDO
	1025	ENDIF
	1026	ENDIF
	1027
	1028	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
	1029
	1030	!
	1031	!-- w-velocity component
	1032	CALL cpu_log( log_point(7), 'w-equation', 'start' )
	1033
[1001]	1034	tend = 0.0
	1035	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	1036	IF ( ws_scheme_mom ) THEN
	1037	CALL advec_w_ws
	1038	ELSE
	1039	CALL advec_w_pw
	1040	ENDIF
	1041	ELSE
[1001]	1042	CALL advec_w_up
[736]	1043	ENDIF
[1001]	1044	CALL diffusion_w
[736]	1045	CALL coriolis( 3 )
[940]	1046
	1047	IF ( .NOT. neutral ) THEN
	1048	IF ( ocean ) THEN
[1179]	1049	CALL buoyancy( rho, 3 )
[736]	1050	ELSE
[940]	1051	IF ( .NOT. humidity ) THEN
[1179]	1052	CALL buoyancy( pt, 3 )
[940]	1053	ELSE
[1179]	1054	CALL buoyancy( vpt, 3 )
[940]	1055	ENDIF
[736]	1056	ENDIF
	1057	ENDIF
	1058
	1059	!
	1060	!-- Drag by plant canopy
	1061	IF ( plant_canopy ) CALL plant_canopy_model( 3 )
	1062
	1063	CALL user_actions( 'w-tendency' )
	1064
	1065	!
	1066	!-- Prognostic equation for w-velocity component
	1067	DO i = nxl, nxr
	1068	DO j = nys, nyn
	1069	DO k = nzb_w_inner(j,i)+1, nzt-1
[1001]	1070	w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	1071	tsc(3) * tw_m(k,j,i) ) &
	1072	- tsc(5) * rdf(k) * w(k,j,i)
[736]	1073	ENDDO
	1074	ENDDO
	1075	ENDDO
	1076
	1077	!
	1078	!-- Calculate tendencies for the next Runge-Kutta step
	1079	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1080	IF ( intermediate_timestep_count == 1 ) THEN
	1081	DO i = nxl, nxr
	1082	DO j = nys, nyn
	1083	DO k = nzb_w_inner(j,i)+1, nzt-1
	1084	tw_m(k,j,i) = tend(k,j,i)
	1085	ENDDO
	1086	ENDDO
	1087	ENDDO
	1088	ELSEIF ( intermediate_timestep_count < &
	1089	intermediate_timestep_count_max ) THEN
	1090	DO i = nxl, nxr
	1091	DO j = nys, nyn
	1092	DO k = nzb_w_inner(j,i)+1, nzt-1
	1093	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
	1094	ENDDO
	1095	ENDDO
	1096	ENDDO
	1097	ENDIF
	1098	ENDIF
	1099
	1100	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
	1101
[940]	1102
[736]	1103	!
[940]	1104	!-- If required, compute prognostic equation for potential temperature
	1105	IF ( .NOT. neutral ) THEN
[736]	1106
[940]	1107	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
	1108
[736]	1109	!
[940]	1110	!-- pt-tendency terms with communication
	1111	sbt = tsc(2)
	1112	IF ( scalar_advec == 'bc-scheme' ) THEN
[736]	1113
[940]	1114	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
[736]	1115	!
[1001]	1116	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[940]	1117	sbt = 1.0
	1118	ENDIF
[736]	1119	tend = 0.0
[940]	1120	CALL advec_s_bc( pt, 'pt' )
[1001]	1121
[736]	1122	ENDIF
[940]	1123
	1124	!
	1125	!-- pt-tendency terms with no communication
[1001]	1126	IF ( scalar_advec /= 'bc-scheme' ) THEN
	1127	tend = 0.0
	1128	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[940]	1129	IF ( ws_scheme_sca ) THEN
	1130	CALL advec_s_ws( pt, 'pt' )
	1131	ELSE
	1132	CALL advec_s_pw( pt )
	1133	ENDIF
	1134	ELSE
[1001]	1135	CALL advec_s_up( pt )
[940]	1136	ENDIF
[736]	1137	ENDIF
	1138
[1001]	1139	CALL diffusion_s( pt, shf, tswst, wall_heatflux )
	1140
[736]	1141	!
[940]	1142	!-- If required compute heating/cooling due to long wave radiation processes
	1143	IF ( radiation ) THEN
	1144	CALL calc_radiation
	1145	ENDIF
[736]	1146
	1147	!
[940]	1148	!-- If required compute impact of latent heat due to precipitation
	1149	IF ( precipitation ) THEN
	1150	CALL impact_of_latent_heat
	1151	ENDIF
[736]	1152
	1153	!
[940]	1154	!-- Consideration of heat sources within the plant canopy
	1155	IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
	1156	CALL plant_canopy_model( 4 )
	1157	ENDIF
[736]	1158
[940]	1159	!
	1160	!-- If required compute influence of large-scale subsidence/ascent
	1161	IF ( large_scale_subsidence ) THEN
	1162	CALL subsidence( tend, pt, pt_init )
	1163	ENDIF
[736]	1164
[1241]	1165	!
	1166	!-- Nudging
	1167	IF ( nudging ) CALL nudge( simulated_time, 'pt' )
	1168
[940]	1169	CALL user_actions( 'pt-tendency' )
[736]	1170
	1171	!
[940]	1172	!-- Prognostic equation for potential temperature
	1173	DO i = nxl, nxr
	1174	DO j = nys, nyn
	1175	DO k = nzb_s_inner(j,i)+1, nzt
[1001]	1176	pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1177	tsc(3) * tpt_m(k,j,i) ) &
	1178	- tsc(5) * ( pt(k,j,i) - pt_init(k) ) *&
	1179	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )
[940]	1180	ENDDO
[736]	1181	ENDDO
	1182	ENDDO
	1183
	1184	!
[940]	1185	!-- Calculate tendencies for the next Runge-Kutta step
	1186	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1187	IF ( intermediate_timestep_count == 1 ) THEN
	1188	DO i = nxl, nxr
	1189	DO j = nys, nyn
	1190	DO k = nzb_s_inner(j,i)+1, nzt
	1191	tpt_m(k,j,i) = tend(k,j,i)
	1192	ENDDO
[736]	1193	ENDDO
	1194	ENDDO
[940]	1195	ELSEIF ( intermediate_timestep_count < &
	1196	intermediate_timestep_count_max ) THEN
	1197	DO i = nxl, nxr
	1198	DO j = nys, nyn
	1199	DO k = nzb_s_inner(j,i)+1, nzt
	1200	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	1201	5.3125 * tpt_m(k,j,i)
	1202	ENDDO
[736]	1203	ENDDO
	1204	ENDDO
[940]	1205	ENDIF
[736]	1206	ENDIF
[940]	1207
	1208	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
	1209
[736]	1210	ENDIF
	1211
	1212	!
	1213	!-- If required, compute prognostic equation for salinity
	1214	IF ( ocean ) THEN
	1215
	1216	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
	1217
	1218	!
	1219	!-- sa-tendency terms with communication
	1220	sbt = tsc(2)
	1221	IF ( scalar_advec == 'bc-scheme' ) THEN
	1222
	1223	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1224	!
[1001]	1225	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[736]	1226	sbt = 1.0
	1227	ENDIF
	1228	tend = 0.0
	1229	CALL advec_s_bc( sa, 'sa' )
[1001]	1230
[736]	1231	ENDIF
	1232
	1233	!
	1234	!-- sa-tendency terms with no communication
[1001]	1235	IF ( scalar_advec /= 'bc-scheme' ) THEN
	1236	tend = 0.0
	1237	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	1238	IF ( ws_scheme_sca ) THEN
	1239	CALL advec_s_ws( sa, 'sa' )
	1240	ELSE
	1241	CALL advec_s_pw( sa )
	1242	ENDIF
	1243	ELSE
[1001]	1244	CALL advec_s_up( sa )
[736]	1245	ENDIF
	1246	ENDIF
[1001]	1247
	1248	CALL diffusion_s( sa, saswsb, saswst, wall_salinityflux )
[736]	1249
	1250	CALL user_actions( 'sa-tendency' )
	1251
	1252	!
	1253	!-- Prognostic equation for salinity
	1254	DO i = nxl, nxr
	1255	DO j = nys, nyn
	1256	DO k = nzb_s_inner(j,i)+1, nzt
[1001]	1257	sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1258	tsc(3) * tsa_m(k,j,i) ) &
	1259	- tsc(5) * rdf_sc(k) * &
	1260	( sa(k,j,i) - sa_init(k) )
[736]	1261	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
	1262	ENDDO
	1263	ENDDO
	1264	ENDDO
	1265
	1266	!
	1267	!-- Calculate tendencies for the next Runge-Kutta step
	1268	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1269	IF ( intermediate_timestep_count == 1 ) THEN
	1270	DO i = nxl, nxr
	1271	DO j = nys, nyn
	1272	DO k = nzb_s_inner(j,i)+1, nzt
	1273	tsa_m(k,j,i) = tend(k,j,i)
	1274	ENDDO
	1275	ENDDO
	1276	ENDDO
	1277	ELSEIF ( intermediate_timestep_count < &
	1278	intermediate_timestep_count_max ) THEN
	1279	DO i = nxl, nxr
	1280	DO j = nys, nyn
	1281	DO k = nzb_s_inner(j,i)+1, nzt
	1282	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	1283	5.3125 * tsa_m(k,j,i)
	1284	ENDDO
	1285	ENDDO
	1286	ENDDO
	1287	ENDIF
	1288	ENDIF
	1289
	1290	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
	1291
	1292	!
	1293	!-- Calculate density by the equation of state for seawater
	1294	CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
	1295	CALL eqn_state_seawater
	1296	CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
	1297
	1298	ENDIF
	1299
	1300	!
	1301	!-- If required, compute prognostic equation for total water content / scalar
	1302	IF ( humidity .OR. passive_scalar ) THEN
	1303
	1304	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
	1305
	1306	!
	1307	!-- Scalar/q-tendency terms with communication
	1308	sbt = tsc(2)
	1309	IF ( scalar_advec == 'bc-scheme' ) THEN
	1310
	1311	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1312	!
[1001]	1313	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[736]	1314	sbt = 1.0
	1315	ENDIF
	1316	tend = 0.0
	1317	CALL advec_s_bc( q, 'q' )
[1001]	1318
[736]	1319	ENDIF
	1320
	1321	!
	1322	!-- Scalar/q-tendency terms with no communication
[1001]	1323	IF ( scalar_advec /= 'bc-scheme' ) THEN
	1324	tend = 0.0
	1325	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	1326	IF ( ws_scheme_sca ) THEN
	1327	CALL advec_s_ws( q, 'q' )
	1328	ELSE
	1329	CALL advec_s_pw( q )
	1330	ENDIF
	1331	ELSE
[1001]	1332	CALL advec_s_up( q )
[736]	1333	ENDIF
	1334	ENDIF
[1001]	1335
	1336	CALL diffusion_s( q, qsws, qswst, wall_qflux )
[736]	1337
	1338	!
	1339	!-- If required compute decrease of total water content due to
	1340	!-- precipitation
	1341	IF ( precipitation ) THEN
	1342	CALL calc_precipitation
	1343	ENDIF
	1344
	1345	!
	1346	!-- Sink or source of scalar concentration due to canopy elements
	1347	IF ( plant_canopy ) CALL plant_canopy_model( 5 )
	1348
	1349	!
	1350	!-- If required compute influence of large-scale subsidence/ascent
[940]	1351	IF ( large_scale_subsidence ) THEN
	1352	CALL subsidence( tend, q, q_init )
[736]	1353	ENDIF
	1354
[1241]	1355	!
	1356	!-- Nudging
	1357	IF ( nudging ) CALL nudge( simulated_time, 'q' )
	1358
[736]	1359	CALL user_actions( 'q-tendency' )
	1360
	1361	!
	1362	!-- Prognostic equation for total water content / scalar
	1363	DO i = nxl, nxr
	1364	DO j = nys, nyn
	1365	DO k = nzb_s_inner(j,i)+1, nzt
[1001]	1366	q_p(k,j,i) = q(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1367	tsc(3) * tq_m(k,j,i) ) &
	1368	- tsc(5) * rdf_sc(k) * &
	1369	( q(k,j,i) - q_init(k) )
[736]	1370	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
	1371	ENDDO
	1372	ENDDO
	1373	ENDDO
	1374
	1375	!
	1376	!-- Calculate tendencies for the next Runge-Kutta step
	1377	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1378	IF ( intermediate_timestep_count == 1 ) THEN
	1379	DO i = nxl, nxr
	1380	DO j = nys, nyn
	1381	DO k = nzb_s_inner(j,i)+1, nzt
	1382	tq_m(k,j,i) = tend(k,j,i)
	1383	ENDDO
	1384	ENDDO
	1385	ENDDO
	1386	ELSEIF ( intermediate_timestep_count < &
	1387	intermediate_timestep_count_max ) THEN
	1388	DO i = nxl, nxr
	1389	DO j = nys, nyn
	1390	DO k = nzb_s_inner(j,i)+1, nzt
	1391	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
	1392	ENDDO
	1393	ENDDO
	1394	ENDDO
	1395	ENDIF
	1396	ENDIF
	1397
	1398	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
	1399
	1400	ENDIF
	1401
	1402	!
	1403	!-- If required, compute prognostic equation for turbulent kinetic
	1404	!-- energy (TKE)
	1405	IF ( .NOT. constant_diffusion ) THEN
	1406
	1407	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
	1408
	1409	sbt = tsc(2)
	1410	IF ( .NOT. use_upstream_for_tke ) THEN
	1411	IF ( scalar_advec == 'bc-scheme' ) THEN
	1412
	1413	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1414	!
[1001]	1415	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[736]	1416	sbt = 1.0
	1417	ENDIF
	1418	tend = 0.0
	1419	CALL advec_s_bc( e, 'e' )
[1001]	1420
[736]	1421	ENDIF
	1422	ENDIF
	1423
	1424	!
	1425	!-- TKE-tendency terms with no communication
[1001]	1426	IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN
[736]	1427	IF ( use_upstream_for_tke ) THEN
	1428	tend = 0.0
	1429	CALL advec_s_up( e )
	1430	ELSE
[1001]	1431	tend = 0.0
	1432	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	1433	IF ( ws_scheme_sca ) THEN
	1434	CALL advec_s_ws( e, 'e' )
	1435	ELSE
	1436	CALL advec_s_pw( e )
	1437	ENDIF
	1438	ELSE
[1001]	1439	CALL advec_s_up( e )
[736]	1440	ENDIF
	1441	ENDIF
[1001]	1442	ENDIF
	1443
	1444	IF ( .NOT. humidity ) THEN
	1445	IF ( ocean ) THEN
	1446	CALL diffusion_e( prho, prho_reference )
[736]	1447	ELSE
[1001]	1448	CALL diffusion_e( pt, pt_reference )
[736]	1449	ENDIF
[1001]	1450	ELSE
	1451	CALL diffusion_e( vpt, pt_reference )
[736]	1452	ENDIF
[1001]	1453
[736]	1454	CALL production_e
	1455
	1456	!
	1457	!-- Additional sink term for flows through plant canopies
	1458	IF ( plant_canopy ) CALL plant_canopy_model( 6 )
	1459	CALL user_actions( 'e-tendency' )
	1460
	1461	!
	1462	!-- Prognostic equation for TKE.
	1463	!-- Eliminate negative TKE values, which can occur due to numerical
	1464	!-- reasons in the course of the integration. In such cases the old TKE
	1465	!-- value is reduced by 90%.
	1466	DO i = nxl, nxr
	1467	DO j = nys, nyn
	1468	DO k = nzb_s_inner(j,i)+1, nzt
[1001]	1469	e_p(k,j,i) = e(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1470	tsc(3) * te_m(k,j,i) )
[736]	1471	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
	1472	ENDDO
	1473	ENDDO
	1474	ENDDO
	1475
	1476	!
	1477	!-- Calculate tendencies for the next Runge-Kutta step
	1478	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1479	IF ( intermediate_timestep_count == 1 ) THEN
	1480	DO i = nxl, nxr
	1481	DO j = nys, nyn
	1482	DO k = nzb_s_inner(j,i)+1, nzt
	1483	te_m(k,j,i) = tend(k,j,i)
	1484	ENDDO
	1485	ENDDO
	1486	ENDDO
	1487	ELSEIF ( intermediate_timestep_count < &
	1488	intermediate_timestep_count_max ) THEN
	1489	DO i = nxl, nxr
	1490	DO j = nys, nyn
	1491	DO k = nzb_s_inner(j,i)+1, nzt
	1492	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
	1493	ENDDO
	1494	ENDDO
	1495	ENDDO
	1496	ENDIF
	1497	ENDIF
	1498
	1499	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
	1500
	1501	ENDIF
	1502
	1503
	1504	END SUBROUTINE prognostic_equations_vector
	1505
	1506
[1015]	1507	SUBROUTINE prognostic_equations_acc
	1508
	1509	!------------------------------------------------------------------------------!
	1510	! Version for accelerator boards
	1511	!------------------------------------------------------------------------------!
	1512
	1513	IMPLICIT NONE
	1514
[1320]	1515	INTEGER(iwp) :: i !:
	1516	INTEGER(iwp) :: j !:
	1517	INTEGER(iwp) :: k !:
	1518	INTEGER(iwp) :: runge_step !:
[1015]	1519
[1320]	1520	REAL(wp) :: sbt !:
	1521
[1015]	1522	!
	1523	!-- Set switch for intermediate Runge-Kutta step
	1524	runge_step = 0
	1525	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1526	IF ( intermediate_timestep_count == 1 ) THEN
	1527	runge_step = 1
	1528	ELSEIF ( intermediate_timestep_count < &
	1529	intermediate_timestep_count_max ) THEN
	1530	runge_step = 2
	1531	ENDIF
	1532	ENDIF
	1533
	1534	!
	1535	!-- u-velocity component
	1536	!++ Statistics still not ported to accelerators
[1179]	1537	!$acc update device( hom, ref_state )
[1015]	1538	CALL cpu_log( log_point(5), 'u-equation', 'start' )
	1539
	1540	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1541	IF ( ws_scheme_mom ) THEN
	1542	CALL advec_u_ws_acc
	1543	ELSE
	1544	tend = 0.0 ! to be removed later??
	1545	CALL advec_u_pw
	1546	ENDIF
	1547	ELSE
	1548	CALL advec_u_up
	1549	ENDIF
	1550	CALL diffusion_u_acc
	1551	CALL coriolis_acc( 1 )
	1552	IF ( sloping_surface .AND. .NOT. neutral ) THEN
[1179]	1553	CALL buoyancy( pt, 1 )
[1015]	1554	ENDIF
	1555
	1556	!
	1557	!-- Drag by plant canopy
	1558	IF ( plant_canopy ) CALL plant_canopy_model( 1 )
	1559
	1560	!
	1561	!-- External pressure gradient
	1562	IF ( dp_external ) THEN
[1128]	1563	DO i = i_left, i_right
	1564	DO j = j_south, j_north
[1015]	1565	DO k = dp_level_ind_b+1, nzt
	1566	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
	1567	ENDDO
	1568	ENDDO
	1569	ENDDO
	1570	ENDIF
	1571
[1246]	1572	!
	1573	!-- Nudging
	1574	IF ( nudging ) CALL nudge( simulated_time, 'u' )
	1575
[1015]	1576	CALL user_actions( 'u-tendency' )
	1577
	1578	!
	1579	!-- Prognostic equation for u-velocity component
	1580	!$acc kernels present( nzb_u_inner, rdf, tend, tu_m, u, ug, u_p )
[1257]	1581	!$acc loop independent
[1128]	1582	DO i = i_left, i_right
[1257]	1583	!$acc loop independent
[1128]	1584	DO j = j_south, j_north
[1257]	1585	!$acc loop independent
[1015]	1586	DO k = 1, nzt
	1587	IF ( k > nzb_u_inner(j,i) ) THEN
	1588	u_p(k,j,i) = u(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	1589	tsc(3) * tu_m(k,j,i) ) &
	1590	- tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
	1591	!
	1592	!-- Tendencies for the next Runge-Kutta step
	1593	IF ( runge_step == 1 ) THEN
	1594	tu_m(k,j,i) = tend(k,j,i)
	1595	ELSEIF ( runge_step == 2 ) THEN
	1596	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
	1597	ENDIF
	1598	ENDIF
	1599	ENDDO
	1600	ENDDO
	1601	ENDDO
	1602	!$acc end kernels
	1603
	1604	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
	1605
	1606	!
	1607	!-- v-velocity component
	1608	CALL cpu_log( log_point(6), 'v-equation', 'start' )
	1609
	1610	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1611	IF ( ws_scheme_mom ) THEN
	1612	CALL advec_v_ws_acc
	1613	ELSE
	1614	tend = 0.0 ! to be removed later??
	1615	CALL advec_v_pw
	1616	END IF
	1617	ELSE
	1618	CALL advec_v_up
	1619	ENDIF
	1620	CALL diffusion_v_acc
	1621	CALL coriolis_acc( 2 )
	1622
	1623	!
	1624	!-- Drag by plant canopy
	1625	IF ( plant_canopy ) CALL plant_canopy_model( 2 )
	1626
	1627	!
	1628	!-- External pressure gradient
	1629	IF ( dp_external ) THEN
[1128]	1630	DO i = i_left, i_right
	1631	DO j = j_south, j_north
[1015]	1632	DO k = dp_level_ind_b+1, nzt
	1633	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
	1634	ENDDO
	1635	ENDDO
	1636	ENDDO
	1637	ENDIF
	1638
[1246]	1639	!
	1640	!-- Nudging
	1641	IF ( nudging ) CALL nudge( simulated_time, 'v' )
	1642
[1015]	1643	CALL user_actions( 'v-tendency' )
	1644
	1645	!
	1646	!-- Prognostic equation for v-velocity component
	1647	!$acc kernels present( nzb_v_inner, rdf, tend, tv_m, v, vg, v_p )
[1257]	1648	!$acc loop independent
[1128]	1649	DO i = i_left, i_right
[1257]	1650	!$acc loop independent
[1128]	1651	DO j = j_south, j_north
[1257]	1652	!$acc loop independent
[1015]	1653	DO k = 1, nzt
	1654	IF ( k > nzb_v_inner(j,i) ) THEN
	1655	v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	1656	tsc(3) * tv_m(k,j,i) ) &
	1657	- tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
	1658	!
	1659	!-- Tendencies for the next Runge-Kutta step
	1660	IF ( runge_step == 1 ) THEN
	1661	tv_m(k,j,i) = tend(k,j,i)
	1662	ELSEIF ( runge_step == 2 ) THEN
	1663	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
	1664	ENDIF
	1665	ENDIF
	1666	ENDDO
	1667	ENDDO
	1668	ENDDO
	1669	!$acc end kernels
	1670
	1671	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
	1672
	1673	!
	1674	!-- w-velocity component
	1675	CALL cpu_log( log_point(7), 'w-equation', 'start' )
	1676
	1677	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1678	IF ( ws_scheme_mom ) THEN
	1679	CALL advec_w_ws_acc
	1680	ELSE
	1681	tend = 0.0 ! to be removed later??
	1682	CALL advec_w_pw
	1683	ENDIF
	1684	ELSE
	1685	CALL advec_w_up
	1686	ENDIF
	1687	CALL diffusion_w_acc
	1688	CALL coriolis_acc( 3 )
	1689
	1690	IF ( .NOT. neutral ) THEN
	1691	IF ( ocean ) THEN
[1179]	1692	CALL buoyancy( rho, 3 )
[1015]	1693	ELSE
	1694	IF ( .NOT. humidity ) THEN
[1179]	1695	CALL buoyancy_acc( pt, 3 )
[1015]	1696	ELSE
[1179]	1697	CALL buoyancy( vpt, 3 )
[1015]	1698	ENDIF
	1699	ENDIF
	1700	ENDIF
	1701
	1702	!
	1703	!-- Drag by plant canopy
	1704	IF ( plant_canopy ) CALL plant_canopy_model( 3 )
	1705
	1706	CALL user_actions( 'w-tendency' )
	1707
	1708	!
	1709	!-- Prognostic equation for w-velocity component
	1710	!$acc kernels present( nzb_w_inner, rdf, tend, tw_m, w, w_p )
[1257]	1711	!$acc loop independent
[1128]	1712	DO i = i_left, i_right
[1257]	1713	!$acc loop independent
[1128]	1714	DO j = j_south, j_north
[1257]	1715	!$acc loop independent
[1015]	1716	DO k = 1, nzt-1
	1717	IF ( k > nzb_w_inner(j,i) ) THEN
	1718	w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	1719	tsc(3) * tw_m(k,j,i) ) &
	1720	- tsc(5) * rdf(k) * w(k,j,i)
	1721	!
	1722	!-- Tendencies for the next Runge-Kutta step
	1723	IF ( runge_step == 1 ) THEN
	1724	tw_m(k,j,i) = tend(k,j,i)
	1725	ELSEIF ( runge_step == 2 ) THEN
	1726	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
	1727	ENDIF
	1728	ENDIF
	1729	ENDDO
	1730	ENDDO
	1731	ENDDO
	1732	!$acc end kernels
	1733
	1734	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
	1735
	1736
	1737	!
	1738	!-- If required, compute prognostic equation for potential temperature
	1739	IF ( .NOT. neutral ) THEN
	1740
	1741	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
	1742
	1743	!
	1744	!-- pt-tendency terms with communication
	1745	sbt = tsc(2)
	1746	IF ( scalar_advec == 'bc-scheme' ) THEN
	1747
	1748	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1749	!
	1750	!-- Bott-Chlond scheme always uses Euler time step. Thus:
	1751	sbt = 1.0
	1752	ENDIF
	1753	tend = 0.0
	1754	CALL advec_s_bc( pt, 'pt' )
	1755
	1756	ENDIF
	1757
	1758	!
	1759	!-- pt-tendency terms with no communication
	1760	IF ( scalar_advec /= 'bc-scheme' ) THEN
	1761	tend = 0.0
	1762	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1763	IF ( ws_scheme_sca ) THEN
	1764	CALL advec_s_ws_acc( pt, 'pt' )
	1765	ELSE
	1766	tend = 0.0 ! to be removed later??
	1767	CALL advec_s_pw( pt )
	1768	ENDIF
	1769	ELSE
	1770	CALL advec_s_up( pt )
	1771	ENDIF
	1772	ENDIF
	1773
	1774	CALL diffusion_s_acc( pt, shf, tswst, wall_heatflux )
	1775
	1776	!
	1777	!-- If required compute heating/cooling due to long wave radiation processes
	1778	IF ( radiation ) THEN
	1779	CALL calc_radiation
	1780	ENDIF
	1781
	1782	!
	1783	!-- If required compute impact of latent heat due to precipitation
	1784	IF ( precipitation ) THEN
	1785	CALL impact_of_latent_heat
	1786	ENDIF
	1787
	1788	!
	1789	!-- Consideration of heat sources within the plant canopy
	1790	IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
	1791	CALL plant_canopy_model( 4 )
	1792	ENDIF
	1793
	1794	!
	1795	!-- If required compute influence of large-scale subsidence/ascent
	1796	IF ( large_scale_subsidence ) THEN
	1797	CALL subsidence( tend, pt, pt_init )
	1798	ENDIF
	1799
[1246]	1800	!
	1801	!-- Nudging
	1802	IF ( nudging ) CALL nudge( simulated_time, 'pt' )
	1803
[1015]	1804	CALL user_actions( 'pt-tendency' )
	1805
	1806	!
	1807	!-- Prognostic equation for potential temperature
	1808	!$acc kernels present( nzb_s_inner, rdf_sc, ptdf_x, ptdf_y, pt_init ) &
	1809	!$acc present( tend, tpt_m, pt, pt_p )
[1257]	1810	!$acc loop independent
[1128]	1811	DO i = i_left, i_right
[1257]	1812	!$acc loop independent
[1128]	1813	DO j = j_south, j_north
[1257]	1814	!$acc loop independent
[1015]	1815	DO k = 1, nzt
	1816	IF ( k > nzb_s_inner(j,i) ) THEN
	1817	pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1818	tsc(3) * tpt_m(k,j,i) ) &
	1819	- tsc(5) * ( pt(k,j,i) - pt_init(k) ) *&
	1820	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )
	1821	!
	1822	!-- Tendencies for the next Runge-Kutta step
	1823	IF ( runge_step == 1 ) THEN
	1824	tpt_m(k,j,i) = tend(k,j,i)
	1825	ELSEIF ( runge_step == 2 ) THEN
	1826	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
	1827	ENDIF
	1828	ENDIF
	1829	ENDDO
	1830	ENDDO
	1831	ENDDO
	1832	!$acc end kernels
	1833
	1834	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
	1835
	1836	ENDIF
	1837
	1838	!
	1839	!-- If required, compute prognostic equation for salinity
	1840	IF ( ocean ) THEN
	1841
	1842	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
	1843
	1844	!
	1845	!-- sa-tendency terms with communication
	1846	sbt = tsc(2)
	1847	IF ( scalar_advec == 'bc-scheme' ) THEN
	1848
	1849	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1850	!
	1851	!-- Bott-Chlond scheme always uses Euler time step. Thus:
	1852	sbt = 1.0
	1853	ENDIF
	1854	tend = 0.0
	1855	CALL advec_s_bc( sa, 'sa' )
	1856
	1857	ENDIF
	1858
	1859	!
	1860	!-- sa-tendency terms with no communication
	1861	IF ( scalar_advec /= 'bc-scheme' ) THEN
	1862	tend = 0.0
	1863	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1864	IF ( ws_scheme_sca ) THEN
	1865	CALL advec_s_ws( sa, 'sa' )
	1866	ELSE
	1867	CALL advec_s_pw( sa )
	1868	ENDIF
	1869	ELSE
	1870	CALL advec_s_up( sa )
	1871	ENDIF
	1872	ENDIF
	1873
	1874	CALL diffusion_s( sa, saswsb, saswst, wall_salinityflux )
	1875
	1876	CALL user_actions( 'sa-tendency' )
	1877
	1878	!
	1879	!-- Prognostic equation for salinity
[1128]	1880	DO i = i_left, i_right
	1881	DO j = j_south, j_north
[1015]	1882	DO k = nzb_s_inner(j,i)+1, nzt
	1883	sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1884	tsc(3) * tsa_m(k,j,i) ) &
	1885	- tsc(5) * rdf_sc(k) * &
	1886	( sa(k,j,i) - sa_init(k) )
	1887	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
	1888	!
	1889	!-- Tendencies for the next Runge-Kutta step
	1890	IF ( runge_step == 1 ) THEN
	1891	tsa_m(k,j,i) = tend(k,j,i)
	1892	ELSEIF ( runge_step == 2 ) THEN
	1893	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tsa_m(k,j,i)
	1894	ENDIF
	1895	ENDDO
	1896	ENDDO
	1897	ENDDO
	1898
	1899	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
	1900
	1901	!
	1902	!-- Calculate density by the equation of state for seawater
	1903	CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
	1904	CALL eqn_state_seawater
	1905	CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
	1906
	1907	ENDIF
	1908
	1909	!
	1910	!-- If required, compute prognostic equation for total water content / scalar
	1911	IF ( humidity .OR. passive_scalar ) THEN
	1912
	1913	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
	1914
	1915	!
	1916	!-- Scalar/q-tendency terms with communication
	1917	sbt = tsc(2)
	1918	IF ( scalar_advec == 'bc-scheme' ) THEN
	1919
	1920	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1921	!
	1922	!-- Bott-Chlond scheme always uses Euler time step. Thus:
	1923	sbt = 1.0
	1924	ENDIF
	1925	tend = 0.0
	1926	CALL advec_s_bc( q, 'q' )
	1927
	1928	ENDIF
	1929
	1930	!
	1931	!-- Scalar/q-tendency terms with no communication
	1932	IF ( scalar_advec /= 'bc-scheme' ) THEN
	1933	tend = 0.0
	1934	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1935	IF ( ws_scheme_sca ) THEN
	1936	CALL advec_s_ws( q, 'q' )
	1937	ELSE
	1938	CALL advec_s_pw( q )
	1939	ENDIF
	1940	ELSE
	1941	CALL advec_s_up( q )
	1942	ENDIF
	1943	ENDIF
	1944
	1945	CALL diffusion_s( q, qsws, qswst, wall_qflux )
	1946
	1947	!
	1948	!-- If required compute decrease of total water content due to
	1949	!-- precipitation
	1950	IF ( precipitation ) THEN
	1951	CALL calc_precipitation
	1952	ENDIF
	1953
	1954	!
	1955	!-- Sink or source of scalar concentration due to canopy elements
	1956	IF ( plant_canopy ) CALL plant_canopy_model( 5 )
	1957
	1958	!
	1959	!-- If required compute influence of large-scale subsidence/ascent
	1960	IF ( large_scale_subsidence ) THEN
	1961	CALL subsidence( tend, q, q_init )
	1962	ENDIF
	1963
[1246]	1964	!
	1965	!-- Nudging
	1966	IF ( nudging ) CALL nudge( simulated_time, 'q' )
	1967
[1015]	1968	CALL user_actions( 'q-tendency' )
	1969
	1970	!
	1971	!-- Prognostic equation for total water content / scalar
[1128]	1972	DO i = i_left, i_right
	1973	DO j = j_south, j_north
[1015]	1974	DO k = nzb_s_inner(j,i)+1, nzt
	1975	q_p(k,j,i) = q(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1976	tsc(3) * tq_m(k,j,i) ) &
	1977	- tsc(5) * rdf_sc(k) * &
	1978	( q(k,j,i) - q_init(k) )
	1979	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
	1980	!
	1981	!-- Tendencies for the next Runge-Kutta step
	1982	IF ( runge_step == 1 ) THEN
	1983	tq_m(k,j,i) = tend(k,j,i)
	1984	ELSEIF ( runge_step == 2 ) THEN
	1985	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
	1986	ENDIF
	1987	ENDDO
	1988	ENDDO
	1989	ENDDO
	1990
	1991	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
	1992
	1993	ENDIF
	1994
	1995	!
	1996	!-- If required, compute prognostic equation for turbulent kinetic
	1997	!-- energy (TKE)
	1998	IF ( .NOT. constant_diffusion ) THEN
	1999
	2000	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
	2001
	2002	sbt = tsc(2)
	2003	IF ( .NOT. use_upstream_for_tke ) THEN
	2004	IF ( scalar_advec == 'bc-scheme' ) THEN
	2005
	2006	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	2007	!
	2008	!-- Bott-Chlond scheme always uses Euler time step. Thus:
	2009	sbt = 1.0
	2010	ENDIF
	2011	tend = 0.0
	2012	CALL advec_s_bc( e, 'e' )
	2013
	2014	ENDIF
	2015	ENDIF
	2016
	2017	!
	2018	!-- TKE-tendency terms with no communication
	2019	IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN
	2020	IF ( use_upstream_for_tke ) THEN
	2021	tend = 0.0
	2022	CALL advec_s_up( e )
	2023	ELSE
	2024	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	2025	IF ( ws_scheme_sca ) THEN
	2026	CALL advec_s_ws_acc( e, 'e' )
	2027	ELSE
	2028	tend = 0.0 ! to be removed later??
	2029	CALL advec_s_pw( e )
	2030	ENDIF
	2031	ELSE
	2032	tend = 0.0 ! to be removed later??
	2033	CALL advec_s_up( e )
	2034	ENDIF
	2035	ENDIF
	2036	ENDIF
	2037
	2038	IF ( .NOT. humidity ) THEN
	2039	IF ( ocean ) THEN
	2040	CALL diffusion_e( prho, prho_reference )
	2041	ELSE
	2042	CALL diffusion_e_acc( pt, pt_reference )
	2043	ENDIF
	2044	ELSE
	2045	CALL diffusion_e( vpt, pt_reference )
	2046	ENDIF
	2047
	2048	CALL production_e_acc
	2049
	2050	!
	2051	!-- Additional sink term for flows through plant canopies
	2052	IF ( plant_canopy ) CALL plant_canopy_model( 6 )
	2053	CALL user_actions( 'e-tendency' )
	2054
	2055	!
	2056	!-- Prognostic equation for TKE.
	2057	!-- Eliminate negative TKE values, which can occur due to numerical
	2058	!-- reasons in the course of the integration. In such cases the old TKE
	2059	!-- value is reduced by 90%.
	2060	!$acc kernels present( e, e_p, nzb_s_inner, tend, te_m )
[1257]	2061	!$acc loop independent
[1128]	2062	DO i = i_left, i_right
[1257]	2063	!$acc loop independent
[1128]	2064	DO j = j_south, j_north
[1257]	2065	!$acc loop independent
[1015]	2066	DO k = 1, nzt
	2067	IF ( k > nzb_s_inner(j,i) ) THEN
	2068	e_p(k,j,i) = e(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	2069	tsc(3) * te_m(k,j,i) )
	2070	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
	2071	!
	2072	!-- Tendencies for the next Runge-Kutta step
	2073	IF ( runge_step == 1 ) THEN
	2074	te_m(k,j,i) = tend(k,j,i)
	2075	ELSEIF ( runge_step == 2 ) THEN
	2076	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
	2077	ENDIF
	2078	ENDIF
	2079	ENDDO
	2080	ENDDO
	2081	ENDDO
	2082	!$acc end kernels
	2083
	2084	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
	2085
	2086	ENDIF
	2087
	2088
	2089	END SUBROUTINE prognostic_equations_acc
	2090
	2091
[736]	2092	END MODULE prognostic_equations_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |