Home

Context Navigation

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

Last change on this file since 1115 was 1115, checked in by hoffmann, 11 years ago
optimization of two-moments cloud physics
Property svn:keywords set to `Id`
File size: 66.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |