Home

Context Navigation

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

Last change on this file since 1016 was 1015, checked in by raasch, 12 years ago
Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added. Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed. mask array is set zero for ghost boundaries
Property svn:keywords set to `Id`
File size: 79.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |