Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |