Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |