Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |