Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |