Home

Context Navigation

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

Last change on this file since 1112 was 1112, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 67.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |