Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |