Home

Context Navigation

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

Last change on this file since 1071 was 1054, checked in by hoffmann, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 67.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |