Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |