Home

Context Navigation

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

Last change on this file since 1342 was 1339, checked in by heinze, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 70.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |