Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |