Home

Context Navigation

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

Last change on this file since 1239 was 1239, checked in by heinze, 11 years ago
routines for nudging and large scale forcing from external file added
Property svn:keywords set to `Id`
File size: 65.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |