Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |