Home

Context Navigation

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

Last change on this file since 1092 was 1092, checked in by raasch, 11 years ago
unused variables remove from several routines
Property svn:keywords set to `Id`
File size: 67.1 KB

Rev	Line
[736]	1	MODULE prognostic_equations_mod
	2
[1036]	3	!--------------------------------------------------------------------------------!
	4	! This file is part of PALM.
	5	!
	6	! PALM is free software: you can redistribute it and/or modify it under the terms
	7	! of the GNU General Public License as published by the Free Software Foundation,
	8	! either version 3 of the License, or (at your option) any later version.
	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
	17	! Copyright 1997-2012 Leibniz University Hannover
	18	!--------------------------------------------------------------------------------!
	19	!
[736]	20	! Current revisions:
[1092]	21	! ------------------
	22	! unused variables removed
[1054]	23	!
	24	! Former revisions:
	25	! -----------------
	26	! $Id: prognostic_equations.f90 1092 2013-02-02 11:24:22Z raasch $
	27	!
	28	! 1053 2012-11-13 17:11:03Z hoffmann
[1053]	29	! implementation of two new prognostic equations for rain drop concentration (nr)
	30	! and rain water content (qr)
[979]	31	!
[1053]	32	! currently, only available for cache loop optimization
[1020]	33	!
[1037]	34	! 1036 2012-10-22 13:43:42Z raasch
	35	! code put under GPL (PALM 3.9)
	36	!
[1020]	37	! 1019 2012-09-28 06:46:45Z raasch
	38	! non-optimized version of prognostic_equations removed
	39	!
[1017]	40	! 1015 2012-09-27 09:23:24Z raasch
	41	! new branch prognostic_equations_acc
	42	! OpenACC statements added + code changes required for GPU optimization
	43	!
[1002]	44	! 1001 2012-09-13 14:08:46Z raasch
	45	! all actions concerning leapfrog- and upstream-spline-scheme removed
	46	!
[979]	47	! 978 2012-08-09 08:28:32Z fricke
[978]	48	! km_damp_x and km_damp_y removed in calls of diffusion_u and diffusion_v
	49	! add ptdf_x, ptdf_y for damping the potential temperature at the inflow
	50	! boundary in case of non-cyclic lateral boundaries
	51	! Bugfix: first thread index changes for WS-scheme at the inflow
[736]	52	!
[941]	53	! 940 2012-07-09 14:31:00Z raasch
	54	! temperature equation can be switched off
	55	!
[786]	56	! 785 2011-11-28 09:47:19Z raasch
	57	! new factor rdf_sc allows separate Rayleigh damping of scalars
	58	!
[737]	59	! 736 2011-08-17 14:13:26Z suehring
	60	! Bugfix: determination of first thread index i for WS-scheme
	61	!
[736]	62	! 709 2011-03-30 09:31:40Z raasch
	63	! formatting adjustments
	64	!
	65	! 673 2011-01-18 16:19:48Z suehring
	66	! Consideration of the pressure gradient (steered by tsc(4)) during the time
	67	! integration removed.
	68	!
	69	! 667 2010-12-23 12:06:00Z suehring/gryschka
	70	! Calls of the advection routines with WS5 added.
	71	! Calls of ws_statistics added to set the statistical arrays to zero after each
	72	! time step.
	73	!
	74	! 531 2010-04-21 06:47:21Z heinze
	75	! add call of subsidence in the equation for humidity / passive scalar
	76	!
	77	! 411 2009-12-11 14:15:58Z heinze
	78	! add call of subsidence in the equation for potential temperature
	79	!
	80	! 388 2009-09-23 09:40:33Z raasch
	81	! prho is used instead of rho in diffusion_e,
	82	! external pressure gradient
	83	!
	84	! 153 2008-03-19 09:41:30Z steinfeld
	85	! add call of plant_canopy_model in the prognostic equation for
	86	! the potential temperature and for the passive scalar
	87	!
	88	! 138 2007-11-28 10:03:58Z letzel
	89	! add call of subroutines that evaluate the canopy drag terms,
	90	! add wall_*flux to parameter list of calls of diffusion_s
	91	!
	92	! 106 2007-08-16 14:30:26Z raasch
	93	! +uswst, vswst as arguments in calls of diffusion_u\|v,
	94	! loops for u and v are starting from index nxlu, nysv, respectively (needed
	95	! for non-cyclic boundary conditions)
	96	!
	97	! 97 2007-06-21 08:23:15Z raasch
	98	! prognostic equation for salinity, density is calculated from equation of
	99	! state for seawater and is used for calculation of buoyancy,
	100	! +eqn_state_seawater_mod
	101	! diffusion_e is called with argument rho in case of ocean runs,
	102	! new argument zw in calls of diffusion_e, new argument pt_/prho_reference
	103	! in calls of buoyancy and diffusion_e, calc_mean_pt_profile renamed
	104	! calc_mean_profile
	105	!
	106	! 75 2007-03-22 09:54:05Z raasch
	107	! checking for negative q and limiting for positive values,
	108	! z0 removed from arguments in calls of diffusion_u/v/w, uxrp, vynp eliminated,
	109	! subroutine names changed to .._noopt, .._cache, and .._vector,
	110	! moisture renamed humidity, Bott-Chlond-scheme can be used in the
	111	! _vector-version
	112	!
	113	! 19 2007-02-23 04:53:48Z raasch
	114	! Calculation of e, q, and pt extended for gridpoint nzt,
	115	! handling of given temperature/humidity/scalar fluxes at top surface
	116	!
	117	! RCS Log replace by Id keyword, revision history cleaned up
	118	!
	119	! Revision 1.21 2006/08/04 15:01:07 raasch
	120	! upstream scheme can be forced to be used for tke (use_upstream_for_tke)
	121	! regardless of the timestep scheme used for the other quantities,
	122	! new argument diss in call of diffusion_e
	123	!
	124	! Revision 1.1 2000/04/13 14:56:27 schroeter
	125	! Initial revision
	126	!
	127	!
	128	! Description:
	129	! ------------
	130	! Solving the prognostic equations.
	131	!------------------------------------------------------------------------------!
	132
	133	USE arrays_3d
	134	USE control_parameters
	135	USE cpulog
	136	USE eqn_state_seawater_mod
	137	USE grid_variables
	138	USE indices
	139	USE interfaces
	140	USE pegrid
	141	USE pointer_interfaces
	142	USE statistics
	143	USE advec_ws
	144	USE advec_s_pw_mod
	145	USE advec_s_up_mod
	146	USE advec_u_pw_mod
	147	USE advec_u_up_mod
	148	USE advec_v_pw_mod
	149	USE advec_v_up_mod
	150	USE advec_w_pw_mod
	151	USE advec_w_up_mod
	152	USE buoyancy_mod
	153	USE calc_precipitation_mod
	154	USE calc_radiation_mod
	155	USE coriolis_mod
	156	USE diffusion_e_mod
	157	USE diffusion_s_mod
	158	USE diffusion_u_mod
	159	USE diffusion_v_mod
	160	USE diffusion_w_mod
	161	USE impact_of_latent_heat_mod
[1053]	162	USE microphysics_mod
[736]	163	USE plant_canopy_model_mod
	164	USE production_e_mod
	165	USE subsidence_mod
	166	USE user_actions_mod
	167
	168
	169	PRIVATE
[1019]	170	PUBLIC prognostic_equations_cache, prognostic_equations_vector, &
	171	prognostic_equations_acc
[736]	172
	173	INTERFACE prognostic_equations_cache
	174	MODULE PROCEDURE prognostic_equations_cache
	175	END INTERFACE prognostic_equations_cache
	176
	177	INTERFACE prognostic_equations_vector
	178	MODULE PROCEDURE prognostic_equations_vector
	179	END INTERFACE prognostic_equations_vector
	180
[1015]	181	INTERFACE prognostic_equations_acc
	182	MODULE PROCEDURE prognostic_equations_acc
	183	END INTERFACE prognostic_equations_acc
[736]	184
[1015]	185
[736]	186	CONTAINS
	187
	188
	189	SUBROUTINE prognostic_equations_cache
	190
	191	!------------------------------------------------------------------------------!
	192	! Version with one optimized loop over all equations. It is only allowed to
	193	! be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
	194	!
	195	! Here the calls of most subroutines are embedded in two DO loops over i and j,
	196	! so communication between CPUs is not allowed (does not make sense) within
	197	! these loops.
	198	!
	199	! (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
	200	!------------------------------------------------------------------------------!
	201
	202	IMPLICIT NONE
	203
	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	INTEGER :: i, j, k
[1001]	842	REAL :: sbt
[736]	843
	844	!
	845	!-- Calculate those variables needed in the tendency terms which need
	846	!-- global communication
[940]	847	IF ( .NOT. neutral ) CALL calc_mean_profile( pt, 4 )
	848	IF ( ocean ) CALL calc_mean_profile( rho, 64 )
	849	IF ( humidity ) CALL calc_mean_profile( vpt, 44 )
[736]	850	IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. &
	851	intermediate_timestep_count == 1 ) CALL ws_statistics
	852
	853	!
	854	!-- u-velocity component
	855	CALL cpu_log( log_point(5), 'u-equation', 'start' )
	856
[1001]	857	tend = 0.0
	858	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	859	IF ( ws_scheme_mom ) THEN
	860	CALL advec_u_ws
	861	ELSE
	862	CALL advec_u_pw
	863	ENDIF
	864	ELSE
[1001]	865	CALL advec_u_up
[736]	866	ENDIF
[1001]	867	CALL diffusion_u
[736]	868	CALL coriolis( 1 )
[940]	869	IF ( sloping_surface .AND. .NOT. neutral ) THEN
	870	CALL buoyancy( pt, pt_reference, 1, 4 )
	871	ENDIF
[736]	872
	873	!
	874	!-- Drag by plant canopy
	875	IF ( plant_canopy ) CALL plant_canopy_model( 1 )
	876
	877	!
	878	!-- External pressure gradient
	879	IF ( dp_external ) THEN
	880	DO i = nxlu, nxr
	881	DO j = nys, nyn
	882	DO k = dp_level_ind_b+1, nzt
	883	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
	884	ENDDO
	885	ENDDO
	886	ENDDO
	887	ENDIF
	888
	889	CALL user_actions( 'u-tendency' )
	890
	891	!
	892	!-- Prognostic equation for u-velocity component
	893	DO i = nxlu, nxr
	894	DO j = nys, nyn
	895	DO k = nzb_u_inner(j,i)+1, nzt
[1001]	896	u_p(k,j,i) = u(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	897	tsc(3) * tu_m(k,j,i) ) &
	898	- tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
[736]	899	ENDDO
	900	ENDDO
	901	ENDDO
	902
	903	!
	904	!-- Calculate tendencies for the next Runge-Kutta step
	905	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	906	IF ( intermediate_timestep_count == 1 ) THEN
	907	DO i = nxlu, nxr
	908	DO j = nys, nyn
	909	DO k = nzb_u_inner(j,i)+1, nzt
	910	tu_m(k,j,i) = tend(k,j,i)
	911	ENDDO
	912	ENDDO
	913	ENDDO
	914	ELSEIF ( intermediate_timestep_count < &
	915	intermediate_timestep_count_max ) THEN
	916	DO i = nxlu, nxr
	917	DO j = nys, nyn
	918	DO k = nzb_u_inner(j,i)+1, nzt
	919	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
	920	ENDDO
	921	ENDDO
	922	ENDDO
	923	ENDIF
	924	ENDIF
	925
	926	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
	927
	928	!
	929	!-- v-velocity component
	930	CALL cpu_log( log_point(6), 'v-equation', 'start' )
	931
[1001]	932	tend = 0.0
	933	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	934	IF ( ws_scheme_mom ) THEN
	935	CALL advec_v_ws
	936	ELSE
	937	CALL advec_v_pw
	938	END IF
	939	ELSE
[1001]	940	CALL advec_v_up
[736]	941	ENDIF
[1001]	942	CALL diffusion_v
[736]	943	CALL coriolis( 2 )
	944
	945	!
	946	!-- Drag by plant canopy
	947	IF ( plant_canopy ) CALL plant_canopy_model( 2 )
	948
	949	!
	950	!-- External pressure gradient
	951	IF ( dp_external ) THEN
	952	DO i = nxl, nxr
	953	DO j = nysv, nyn
	954	DO k = dp_level_ind_b+1, nzt
	955	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
	956	ENDDO
	957	ENDDO
	958	ENDDO
	959	ENDIF
	960
	961	CALL user_actions( 'v-tendency' )
	962
	963	!
	964	!-- Prognostic equation for v-velocity component
	965	DO i = nxl, nxr
	966	DO j = nysv, nyn
	967	DO k = nzb_v_inner(j,i)+1, nzt
[1001]	968	v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	969	tsc(3) * tv_m(k,j,i) ) &
	970	- tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
[736]	971	ENDDO
	972	ENDDO
	973	ENDDO
	974
	975	!
	976	!-- Calculate tendencies for the next Runge-Kutta step
	977	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	978	IF ( intermediate_timestep_count == 1 ) THEN
	979	DO i = nxl, nxr
	980	DO j = nysv, nyn
	981	DO k = nzb_v_inner(j,i)+1, nzt
	982	tv_m(k,j,i) = tend(k,j,i)
	983	ENDDO
	984	ENDDO
	985	ENDDO
	986	ELSEIF ( intermediate_timestep_count < &
	987	intermediate_timestep_count_max ) THEN
	988	DO i = nxl, nxr
	989	DO j = nysv, nyn
	990	DO k = nzb_v_inner(j,i)+1, nzt
	991	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
	992	ENDDO
	993	ENDDO
	994	ENDDO
	995	ENDIF
	996	ENDIF
	997
	998	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
	999
	1000	!
	1001	!-- w-velocity component
	1002	CALL cpu_log( log_point(7), 'w-equation', 'start' )
	1003
[1001]	1004	tend = 0.0
	1005	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	1006	IF ( ws_scheme_mom ) THEN
	1007	CALL advec_w_ws
	1008	ELSE
	1009	CALL advec_w_pw
	1010	ENDIF
	1011	ELSE
[1001]	1012	CALL advec_w_up
[736]	1013	ENDIF
[1001]	1014	CALL diffusion_w
[736]	1015	CALL coriolis( 3 )
[940]	1016
	1017	IF ( .NOT. neutral ) THEN
	1018	IF ( ocean ) THEN
	1019	CALL buoyancy( rho, rho_reference, 3, 64 )
[736]	1020	ELSE
[940]	1021	IF ( .NOT. humidity ) THEN
	1022	CALL buoyancy( pt, pt_reference, 3, 4 )
	1023	ELSE
	1024	CALL buoyancy( vpt, pt_reference, 3, 44 )
	1025	ENDIF
[736]	1026	ENDIF
	1027	ENDIF
	1028
	1029	!
	1030	!-- Drag by plant canopy
	1031	IF ( plant_canopy ) CALL plant_canopy_model( 3 )
	1032
	1033	CALL user_actions( 'w-tendency' )
	1034
	1035	!
	1036	!-- Prognostic equation for w-velocity component
	1037	DO i = nxl, nxr
	1038	DO j = nys, nyn
	1039	DO k = nzb_w_inner(j,i)+1, nzt-1
[1001]	1040	w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	1041	tsc(3) * tw_m(k,j,i) ) &
	1042	- tsc(5) * rdf(k) * w(k,j,i)
[736]	1043	ENDDO
	1044	ENDDO
	1045	ENDDO
	1046
	1047	!
	1048	!-- Calculate tendencies for the next Runge-Kutta step
	1049	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1050	IF ( intermediate_timestep_count == 1 ) THEN
	1051	DO i = nxl, nxr
	1052	DO j = nys, nyn
	1053	DO k = nzb_w_inner(j,i)+1, nzt-1
	1054	tw_m(k,j,i) = tend(k,j,i)
	1055	ENDDO
	1056	ENDDO
	1057	ENDDO
	1058	ELSEIF ( intermediate_timestep_count < &
	1059	intermediate_timestep_count_max ) THEN
	1060	DO i = nxl, nxr
	1061	DO j = nys, nyn
	1062	DO k = nzb_w_inner(j,i)+1, nzt-1
	1063	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
	1064	ENDDO
	1065	ENDDO
	1066	ENDDO
	1067	ENDIF
	1068	ENDIF
	1069
	1070	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
	1071
[940]	1072
[736]	1073	!
[940]	1074	!-- If required, compute prognostic equation for potential temperature
	1075	IF ( .NOT. neutral ) THEN
[736]	1076
[940]	1077	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
	1078
[736]	1079	!
[940]	1080	!-- pt-tendency terms with communication
	1081	sbt = tsc(2)
	1082	IF ( scalar_advec == 'bc-scheme' ) THEN
[736]	1083
[940]	1084	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
[736]	1085	!
[1001]	1086	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[940]	1087	sbt = 1.0
	1088	ENDIF
[736]	1089	tend = 0.0
[940]	1090	CALL advec_s_bc( pt, 'pt' )
[1001]	1091
[736]	1092	ENDIF
[940]	1093
	1094	!
	1095	!-- pt-tendency terms with no communication
[1001]	1096	IF ( scalar_advec /= 'bc-scheme' ) THEN
	1097	tend = 0.0
	1098	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[940]	1099	IF ( ws_scheme_sca ) THEN
	1100	CALL advec_s_ws( pt, 'pt' )
	1101	ELSE
	1102	CALL advec_s_pw( pt )
	1103	ENDIF
	1104	ELSE
[1001]	1105	CALL advec_s_up( pt )
[940]	1106	ENDIF
[736]	1107	ENDIF
	1108
[1001]	1109	CALL diffusion_s( pt, shf, tswst, wall_heatflux )
	1110
[736]	1111	!
[940]	1112	!-- If required compute heating/cooling due to long wave radiation processes
	1113	IF ( radiation ) THEN
	1114	CALL calc_radiation
	1115	ENDIF
[736]	1116
	1117	!
[940]	1118	!-- If required compute impact of latent heat due to precipitation
	1119	IF ( precipitation ) THEN
	1120	CALL impact_of_latent_heat
	1121	ENDIF
[736]	1122
	1123	!
[940]	1124	!-- Consideration of heat sources within the plant canopy
	1125	IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
	1126	CALL plant_canopy_model( 4 )
	1127	ENDIF
[736]	1128
[940]	1129	!
	1130	!-- If required compute influence of large-scale subsidence/ascent
	1131	IF ( large_scale_subsidence ) THEN
	1132	CALL subsidence( tend, pt, pt_init )
	1133	ENDIF
[736]	1134
[940]	1135	CALL user_actions( 'pt-tendency' )
[736]	1136
	1137	!
[940]	1138	!-- Prognostic equation for potential temperature
	1139	DO i = nxl, nxr
	1140	DO j = nys, nyn
	1141	DO k = nzb_s_inner(j,i)+1, nzt
[1001]	1142	pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1143	tsc(3) * tpt_m(k,j,i) ) &
	1144	- tsc(5) * ( pt(k,j,i) - pt_init(k) ) *&
	1145	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )
[940]	1146	ENDDO
[736]	1147	ENDDO
	1148	ENDDO
	1149
	1150	!
[940]	1151	!-- Calculate tendencies for the next Runge-Kutta step
	1152	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1153	IF ( intermediate_timestep_count == 1 ) THEN
	1154	DO i = nxl, nxr
	1155	DO j = nys, nyn
	1156	DO k = nzb_s_inner(j,i)+1, nzt
	1157	tpt_m(k,j,i) = tend(k,j,i)
	1158	ENDDO
[736]	1159	ENDDO
	1160	ENDDO
[940]	1161	ELSEIF ( intermediate_timestep_count < &
	1162	intermediate_timestep_count_max ) THEN
	1163	DO i = nxl, nxr
	1164	DO j = nys, nyn
	1165	DO k = nzb_s_inner(j,i)+1, nzt
	1166	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	1167	5.3125 * tpt_m(k,j,i)
	1168	ENDDO
[736]	1169	ENDDO
	1170	ENDDO
[940]	1171	ENDIF
[736]	1172	ENDIF
[940]	1173
	1174	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
	1175
[736]	1176	ENDIF
	1177
	1178	!
	1179	!-- If required, compute prognostic equation for salinity
	1180	IF ( ocean ) THEN
	1181
	1182	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
	1183
	1184	!
	1185	!-- sa-tendency terms with communication
	1186	sbt = tsc(2)
	1187	IF ( scalar_advec == 'bc-scheme' ) THEN
	1188
	1189	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1190	!
[1001]	1191	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[736]	1192	sbt = 1.0
	1193	ENDIF
	1194	tend = 0.0
	1195	CALL advec_s_bc( sa, 'sa' )
[1001]	1196
[736]	1197	ENDIF
	1198
	1199	!
	1200	!-- sa-tendency terms with no communication
[1001]	1201	IF ( scalar_advec /= 'bc-scheme' ) THEN
	1202	tend = 0.0
	1203	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	1204	IF ( ws_scheme_sca ) THEN
	1205	CALL advec_s_ws( sa, 'sa' )
	1206	ELSE
	1207	CALL advec_s_pw( sa )
	1208	ENDIF
	1209	ELSE
[1001]	1210	CALL advec_s_up( sa )
[736]	1211	ENDIF
	1212	ENDIF
[1001]	1213
	1214	CALL diffusion_s( sa, saswsb, saswst, wall_salinityflux )
[736]	1215
	1216	CALL user_actions( 'sa-tendency' )
	1217
	1218	!
	1219	!-- Prognostic equation for salinity
	1220	DO i = nxl, nxr
	1221	DO j = nys, nyn
	1222	DO k = nzb_s_inner(j,i)+1, nzt
[1001]	1223	sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1224	tsc(3) * tsa_m(k,j,i) ) &
	1225	- tsc(5) * rdf_sc(k) * &
	1226	( sa(k,j,i) - sa_init(k) )
[736]	1227	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
	1228	ENDDO
	1229	ENDDO
	1230	ENDDO
	1231
	1232	!
	1233	!-- Calculate tendencies for the next Runge-Kutta step
	1234	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1235	IF ( intermediate_timestep_count == 1 ) THEN
	1236	DO i = nxl, nxr
	1237	DO j = nys, nyn
	1238	DO k = nzb_s_inner(j,i)+1, nzt
	1239	tsa_m(k,j,i) = tend(k,j,i)
	1240	ENDDO
	1241	ENDDO
	1242	ENDDO
	1243	ELSEIF ( intermediate_timestep_count < &
	1244	intermediate_timestep_count_max ) THEN
	1245	DO i = nxl, nxr
	1246	DO j = nys, nyn
	1247	DO k = nzb_s_inner(j,i)+1, nzt
	1248	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
	1249	5.3125 * tsa_m(k,j,i)
	1250	ENDDO
	1251	ENDDO
	1252	ENDDO
	1253	ENDIF
	1254	ENDIF
	1255
	1256	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
	1257
	1258	!
	1259	!-- Calculate density by the equation of state for seawater
	1260	CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
	1261	CALL eqn_state_seawater
	1262	CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
	1263
	1264	ENDIF
	1265
	1266	!
	1267	!-- If required, compute prognostic equation for total water content / scalar
	1268	IF ( humidity .OR. passive_scalar ) THEN
	1269
	1270	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
	1271
	1272	!
	1273	!-- Scalar/q-tendency terms with communication
	1274	sbt = tsc(2)
	1275	IF ( scalar_advec == 'bc-scheme' ) THEN
	1276
	1277	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1278	!
[1001]	1279	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[736]	1280	sbt = 1.0
	1281	ENDIF
	1282	tend = 0.0
	1283	CALL advec_s_bc( q, 'q' )
[1001]	1284
[736]	1285	ENDIF
	1286
	1287	!
	1288	!-- Scalar/q-tendency terms with no communication
[1001]	1289	IF ( scalar_advec /= 'bc-scheme' ) THEN
	1290	tend = 0.0
	1291	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	1292	IF ( ws_scheme_sca ) THEN
	1293	CALL advec_s_ws( q, 'q' )
	1294	ELSE
	1295	CALL advec_s_pw( q )
	1296	ENDIF
	1297	ELSE
[1001]	1298	CALL advec_s_up( q )
[736]	1299	ENDIF
	1300	ENDIF
[1001]	1301
	1302	CALL diffusion_s( q, qsws, qswst, wall_qflux )
[736]	1303
	1304	!
	1305	!-- If required compute decrease of total water content due to
	1306	!-- precipitation
	1307	IF ( precipitation ) THEN
	1308	CALL calc_precipitation
	1309	ENDIF
	1310
	1311	!
	1312	!-- Sink or source of scalar concentration due to canopy elements
	1313	IF ( plant_canopy ) CALL plant_canopy_model( 5 )
	1314
	1315	!
	1316	!-- If required compute influence of large-scale subsidence/ascent
[940]	1317	IF ( large_scale_subsidence ) THEN
	1318	CALL subsidence( tend, q, q_init )
[736]	1319	ENDIF
	1320
	1321	CALL user_actions( 'q-tendency' )
	1322
	1323	!
	1324	!-- Prognostic equation for total water content / scalar
	1325	DO i = nxl, nxr
	1326	DO j = nys, nyn
	1327	DO k = nzb_s_inner(j,i)+1, nzt
[1001]	1328	q_p(k,j,i) = q(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1329	tsc(3) * tq_m(k,j,i) ) &
	1330	- tsc(5) * rdf_sc(k) * &
	1331	( q(k,j,i) - q_init(k) )
[736]	1332	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
	1333	ENDDO
	1334	ENDDO
	1335	ENDDO
	1336
	1337	!
	1338	!-- Calculate tendencies for the next Runge-Kutta step
	1339	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1340	IF ( intermediate_timestep_count == 1 ) THEN
	1341	DO i = nxl, nxr
	1342	DO j = nys, nyn
	1343	DO k = nzb_s_inner(j,i)+1, nzt
	1344	tq_m(k,j,i) = tend(k,j,i)
	1345	ENDDO
	1346	ENDDO
	1347	ENDDO
	1348	ELSEIF ( intermediate_timestep_count < &
	1349	intermediate_timestep_count_max ) THEN
	1350	DO i = nxl, nxr
	1351	DO j = nys, nyn
	1352	DO k = nzb_s_inner(j,i)+1, nzt
	1353	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
	1354	ENDDO
	1355	ENDDO
	1356	ENDDO
	1357	ENDIF
	1358	ENDIF
	1359
	1360	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
	1361
	1362	ENDIF
	1363
	1364	!
	1365	!-- If required, compute prognostic equation for turbulent kinetic
	1366	!-- energy (TKE)
	1367	IF ( .NOT. constant_diffusion ) THEN
	1368
	1369	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
	1370
	1371	!
	1372	!-- TKE-tendency terms with communication
	1373	CALL production_e_init
	1374
	1375	sbt = tsc(2)
	1376	IF ( .NOT. use_upstream_for_tke ) THEN
	1377	IF ( scalar_advec == 'bc-scheme' ) THEN
	1378
	1379	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1380	!
[1001]	1381	!-- Bott-Chlond scheme always uses Euler time step. Thus:
[736]	1382	sbt = 1.0
	1383	ENDIF
	1384	tend = 0.0
	1385	CALL advec_s_bc( e, 'e' )
[1001]	1386
[736]	1387	ENDIF
	1388	ENDIF
	1389
	1390	!
	1391	!-- TKE-tendency terms with no communication
[1001]	1392	IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN
[736]	1393	IF ( use_upstream_for_tke ) THEN
	1394	tend = 0.0
	1395	CALL advec_s_up( e )
	1396	ELSE
[1001]	1397	tend = 0.0
	1398	IF ( timestep_scheme(1:5) == 'runge' ) THEN
[736]	1399	IF ( ws_scheme_sca ) THEN
	1400	CALL advec_s_ws( e, 'e' )
	1401	ELSE
	1402	CALL advec_s_pw( e )
	1403	ENDIF
	1404	ELSE
[1001]	1405	CALL advec_s_up( e )
[736]	1406	ENDIF
	1407	ENDIF
[1001]	1408	ENDIF
	1409
	1410	IF ( .NOT. humidity ) THEN
	1411	IF ( ocean ) THEN
	1412	CALL diffusion_e( prho, prho_reference )
[736]	1413	ELSE
[1001]	1414	CALL diffusion_e( pt, pt_reference )
[736]	1415	ENDIF
[1001]	1416	ELSE
	1417	CALL diffusion_e( vpt, pt_reference )
[736]	1418	ENDIF
[1001]	1419
[736]	1420	CALL production_e
	1421
	1422	!
	1423	!-- Additional sink term for flows through plant canopies
	1424	IF ( plant_canopy ) CALL plant_canopy_model( 6 )
	1425	CALL user_actions( 'e-tendency' )
	1426
	1427	!
	1428	!-- Prognostic equation for TKE.
	1429	!-- Eliminate negative TKE values, which can occur due to numerical
	1430	!-- reasons in the course of the integration. In such cases the old TKE
	1431	!-- value is reduced by 90%.
	1432	DO i = nxl, nxr
	1433	DO j = nys, nyn
	1434	DO k = nzb_s_inner(j,i)+1, nzt
[1001]	1435	e_p(k,j,i) = e(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1436	tsc(3) * te_m(k,j,i) )
[736]	1437	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
	1438	ENDDO
	1439	ENDDO
	1440	ENDDO
	1441
	1442	!
	1443	!-- Calculate tendencies for the next Runge-Kutta step
	1444	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1445	IF ( intermediate_timestep_count == 1 ) THEN
	1446	DO i = nxl, nxr
	1447	DO j = nys, nyn
	1448	DO k = nzb_s_inner(j,i)+1, nzt
	1449	te_m(k,j,i) = tend(k,j,i)
	1450	ENDDO
	1451	ENDDO
	1452	ENDDO
	1453	ELSEIF ( intermediate_timestep_count < &
	1454	intermediate_timestep_count_max ) THEN
	1455	DO i = nxl, nxr
	1456	DO j = nys, nyn
	1457	DO k = nzb_s_inner(j,i)+1, nzt
	1458	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
	1459	ENDDO
	1460	ENDDO
	1461	ENDDO
	1462	ENDIF
	1463	ENDIF
	1464
	1465	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
	1466
	1467	ENDIF
	1468
	1469
	1470	END SUBROUTINE prognostic_equations_vector
	1471
	1472
[1015]	1473	SUBROUTINE prognostic_equations_acc
	1474
	1475	!------------------------------------------------------------------------------!
	1476	! Version for accelerator boards
	1477	!------------------------------------------------------------------------------!
	1478
	1479	IMPLICIT NONE
	1480
	1481	INTEGER :: i, j, k, runge_step
	1482	REAL :: sbt
	1483
	1484	!
	1485	!-- Set switch for intermediate Runge-Kutta step
	1486	runge_step = 0
	1487	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1488	IF ( intermediate_timestep_count == 1 ) THEN
	1489	runge_step = 1
	1490	ELSEIF ( intermediate_timestep_count < &
	1491	intermediate_timestep_count_max ) THEN
	1492	runge_step = 2
	1493	ENDIF
	1494	ENDIF
	1495
	1496	!
	1497	!-- Calculate those variables needed in the tendency terms which need
	1498	!-- global communication
	1499	IF ( .NOT. neutral ) CALL calc_mean_profile( pt, 4 )
	1500	IF ( ocean ) CALL calc_mean_profile( rho, 64 )
	1501	IF ( humidity ) CALL calc_mean_profile( vpt, 44 )
	1502	IF ( ( ws_scheme_mom .OR. ws_scheme_sca ) .AND. &
	1503	intermediate_timestep_count == 1 ) CALL ws_statistics
	1504
	1505	!
	1506	!-- u-velocity component
	1507	!++ Statistics still not ported to accelerators
	1508	!$acc update device( hom )
	1509	CALL cpu_log( log_point(5), 'u-equation', 'start' )
	1510
	1511	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1512	IF ( ws_scheme_mom ) THEN
	1513	CALL advec_u_ws_acc
	1514	ELSE
	1515	tend = 0.0 ! to be removed later??
	1516	CALL advec_u_pw
	1517	ENDIF
	1518	ELSE
	1519	CALL advec_u_up
	1520	ENDIF
	1521	CALL diffusion_u_acc
	1522	CALL coriolis_acc( 1 )
	1523	IF ( sloping_surface .AND. .NOT. neutral ) THEN
	1524	CALL buoyancy( pt, pt_reference, 1, 4 )
	1525	ENDIF
	1526
	1527	!
	1528	!-- Drag by plant canopy
	1529	IF ( plant_canopy ) CALL plant_canopy_model( 1 )
	1530
	1531	!
	1532	!-- External pressure gradient
	1533	IF ( dp_external ) THEN
	1534	DO i = nxlu, nxr
	1535	DO j = nys, nyn
	1536	DO k = dp_level_ind_b+1, nzt
	1537	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
	1538	ENDDO
	1539	ENDDO
	1540	ENDDO
	1541	ENDIF
	1542
	1543	CALL user_actions( 'u-tendency' )
	1544
	1545	!
	1546	!-- Prognostic equation for u-velocity component
	1547	!$acc kernels present( nzb_u_inner, rdf, tend, tu_m, u, ug, u_p )
	1548	!$acc loop
	1549	DO i = nxlu, nxr
	1550	DO j = nys, nyn
	1551	!$acc loop vector( 32 )
	1552	DO k = 1, nzt
	1553	IF ( k > nzb_u_inner(j,i) ) THEN
	1554	u_p(k,j,i) = u(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	1555	tsc(3) * tu_m(k,j,i) ) &
	1556	- tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
	1557	!
	1558	!-- Tendencies for the next Runge-Kutta step
	1559	IF ( runge_step == 1 ) THEN
	1560	tu_m(k,j,i) = tend(k,j,i)
	1561	ELSEIF ( runge_step == 2 ) THEN
	1562	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
	1563	ENDIF
	1564	ENDIF
	1565	ENDDO
	1566	ENDDO
	1567	ENDDO
	1568	!$acc end kernels
	1569
	1570	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
	1571	!$acc update host( u_p )
	1572
	1573	!
	1574	!-- v-velocity component
	1575	CALL cpu_log( log_point(6), 'v-equation', 'start' )
	1576
	1577	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1578	IF ( ws_scheme_mom ) THEN
	1579	CALL advec_v_ws_acc
	1580	ELSE
	1581	tend = 0.0 ! to be removed later??
	1582	CALL advec_v_pw
	1583	END IF
	1584	ELSE
	1585	CALL advec_v_up
	1586	ENDIF
	1587	CALL diffusion_v_acc
	1588	CALL coriolis_acc( 2 )
	1589
	1590	!
	1591	!-- Drag by plant canopy
	1592	IF ( plant_canopy ) CALL plant_canopy_model( 2 )
	1593
	1594	!
	1595	!-- External pressure gradient
	1596	IF ( dp_external ) THEN
	1597	DO i = nxl, nxr
	1598	DO j = nysv, nyn
	1599	DO k = dp_level_ind_b+1, nzt
	1600	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
	1601	ENDDO
	1602	ENDDO
	1603	ENDDO
	1604	ENDIF
	1605
	1606	CALL user_actions( 'v-tendency' )
	1607
	1608	!
	1609	!-- Prognostic equation for v-velocity component
	1610	!$acc kernels present( nzb_v_inner, rdf, tend, tv_m, v, vg, v_p )
	1611	!$acc loop
	1612	DO i = nxl, nxr
	1613	DO j = nysv, nyn
	1614	!$acc loop vector( 32 )
	1615	DO k = 1, nzt
	1616	IF ( k > nzb_v_inner(j,i) ) THEN
	1617	v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	1618	tsc(3) * tv_m(k,j,i) ) &
	1619	- tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
	1620	!
	1621	!-- Tendencies for the next Runge-Kutta step
	1622	IF ( runge_step == 1 ) THEN
	1623	tv_m(k,j,i) = tend(k,j,i)
	1624	ELSEIF ( runge_step == 2 ) THEN
	1625	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
	1626	ENDIF
	1627	ENDIF
	1628	ENDDO
	1629	ENDDO
	1630	ENDDO
	1631	!$acc end kernels
	1632
	1633	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
	1634	!$acc update host( v_p )
	1635
	1636	!
	1637	!-- w-velocity component
	1638	CALL cpu_log( log_point(7), 'w-equation', 'start' )
	1639
	1640	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1641	IF ( ws_scheme_mom ) THEN
	1642	CALL advec_w_ws_acc
	1643	ELSE
	1644	tend = 0.0 ! to be removed later??
	1645	CALL advec_w_pw
	1646	ENDIF
	1647	ELSE
	1648	CALL advec_w_up
	1649	ENDIF
	1650	CALL diffusion_w_acc
	1651	CALL coriolis_acc( 3 )
	1652
	1653	IF ( .NOT. neutral ) THEN
	1654	IF ( ocean ) THEN
	1655	CALL buoyancy( rho, rho_reference, 3, 64 )
	1656	ELSE
	1657	IF ( .NOT. humidity ) THEN
	1658	CALL buoyancy_acc( pt, pt_reference, 3, 4 )
	1659	ELSE
	1660	CALL buoyancy( vpt, pt_reference, 3, 44 )
	1661	ENDIF
	1662	ENDIF
	1663	ENDIF
	1664
	1665	!
	1666	!-- Drag by plant canopy
	1667	IF ( plant_canopy ) CALL plant_canopy_model( 3 )
	1668
	1669	CALL user_actions( 'w-tendency' )
	1670
	1671	!
	1672	!-- Prognostic equation for w-velocity component
	1673	!$acc kernels present( nzb_w_inner, rdf, tend, tw_m, w, w_p )
	1674	!$acc loop
	1675	DO i = nxl, nxr
	1676	DO j = nys, nyn
	1677	!$acc loop vector( 32 )
	1678	DO k = 1, nzt-1
	1679	IF ( k > nzb_w_inner(j,i) ) THEN
	1680	w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
	1681	tsc(3) * tw_m(k,j,i) ) &
	1682	- tsc(5) * rdf(k) * w(k,j,i)
	1683	!
	1684	!-- Tendencies for the next Runge-Kutta step
	1685	IF ( runge_step == 1 ) THEN
	1686	tw_m(k,j,i) = tend(k,j,i)
	1687	ELSEIF ( runge_step == 2 ) THEN
	1688	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
	1689	ENDIF
	1690	ENDIF
	1691	ENDDO
	1692	ENDDO
	1693	ENDDO
	1694	!$acc end kernels
	1695
	1696	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
	1697	!$acc update host( w_p )
	1698
	1699
	1700	!
	1701	!-- If required, compute prognostic equation for potential temperature
	1702	IF ( .NOT. neutral ) THEN
	1703
	1704	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
	1705
	1706	!
	1707	!-- pt-tendency terms with communication
	1708	sbt = tsc(2)
	1709	IF ( scalar_advec == 'bc-scheme' ) THEN
	1710
	1711	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1712	!
	1713	!-- Bott-Chlond scheme always uses Euler time step. Thus:
	1714	sbt = 1.0
	1715	ENDIF
	1716	tend = 0.0
	1717	CALL advec_s_bc( pt, 'pt' )
	1718
	1719	ENDIF
	1720
	1721	!
	1722	!-- pt-tendency terms with no communication
	1723	IF ( scalar_advec /= 'bc-scheme' ) THEN
	1724	tend = 0.0
	1725	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1726	IF ( ws_scheme_sca ) THEN
	1727	CALL advec_s_ws_acc( pt, 'pt' )
	1728	ELSE
	1729	tend = 0.0 ! to be removed later??
	1730	CALL advec_s_pw( pt )
	1731	ENDIF
	1732	ELSE
	1733	CALL advec_s_up( pt )
	1734	ENDIF
	1735	ENDIF
	1736
	1737	CALL diffusion_s_acc( pt, shf, tswst, wall_heatflux )
	1738
	1739	!
	1740	!-- If required compute heating/cooling due to long wave radiation processes
	1741	IF ( radiation ) THEN
	1742	CALL calc_radiation
	1743	ENDIF
	1744
	1745	!
	1746	!-- If required compute impact of latent heat due to precipitation
	1747	IF ( precipitation ) THEN
	1748	CALL impact_of_latent_heat
	1749	ENDIF
	1750
	1751	!
	1752	!-- Consideration of heat sources within the plant canopy
	1753	IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
	1754	CALL plant_canopy_model( 4 )
	1755	ENDIF
	1756
	1757	!
	1758	!-- If required compute influence of large-scale subsidence/ascent
	1759	IF ( large_scale_subsidence ) THEN
	1760	CALL subsidence( tend, pt, pt_init )
	1761	ENDIF
	1762
	1763	CALL user_actions( 'pt-tendency' )
	1764
	1765	!
	1766	!-- Prognostic equation for potential temperature
	1767	!$acc kernels present( nzb_s_inner, rdf_sc, ptdf_x, ptdf_y, pt_init ) &
	1768	!$acc present( tend, tpt_m, pt, pt_p )
	1769	!$acc loop
	1770	DO i = nxl, nxr
	1771	DO j = nys, nyn
	1772	!$acc loop vector( 32 )
	1773	DO k = 1, nzt
	1774	IF ( k > nzb_s_inner(j,i) ) THEN
	1775	pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1776	tsc(3) * tpt_m(k,j,i) ) &
	1777	- tsc(5) * ( pt(k,j,i) - pt_init(k) ) *&
	1778	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )
	1779	!
	1780	!-- Tendencies for the next Runge-Kutta step
	1781	IF ( runge_step == 1 ) THEN
	1782	tpt_m(k,j,i) = tend(k,j,i)
	1783	ELSEIF ( runge_step == 2 ) THEN
	1784	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
	1785	ENDIF
	1786	ENDIF
	1787	ENDDO
	1788	ENDDO
	1789	ENDDO
	1790	!$acc end kernels
	1791
	1792	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
	1793	!$acc update host( pt_p )
	1794
	1795	ENDIF
	1796
	1797	!
	1798	!-- If required, compute prognostic equation for salinity
	1799	IF ( ocean ) THEN
	1800
	1801	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
	1802
	1803	!
	1804	!-- sa-tendency terms with communication
	1805	sbt = tsc(2)
	1806	IF ( scalar_advec == 'bc-scheme' ) THEN
	1807
	1808	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1809	!
	1810	!-- Bott-Chlond scheme always uses Euler time step. Thus:
	1811	sbt = 1.0
	1812	ENDIF
	1813	tend = 0.0
	1814	CALL advec_s_bc( sa, 'sa' )
	1815
	1816	ENDIF
	1817
	1818	!
	1819	!-- sa-tendency terms with no communication
	1820	IF ( scalar_advec /= 'bc-scheme' ) THEN
	1821	tend = 0.0
	1822	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1823	IF ( ws_scheme_sca ) THEN
	1824	CALL advec_s_ws( sa, 'sa' )
	1825	ELSE
	1826	CALL advec_s_pw( sa )
	1827	ENDIF
	1828	ELSE
	1829	CALL advec_s_up( sa )
	1830	ENDIF
	1831	ENDIF
	1832
	1833	CALL diffusion_s( sa, saswsb, saswst, wall_salinityflux )
	1834
	1835	CALL user_actions( 'sa-tendency' )
	1836
	1837	!
	1838	!-- Prognostic equation for salinity
	1839	DO i = nxl, nxr
	1840	DO j = nys, nyn
	1841	DO k = nzb_s_inner(j,i)+1, nzt
	1842	sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1843	tsc(3) * tsa_m(k,j,i) ) &
	1844	- tsc(5) * rdf_sc(k) * &
	1845	( sa(k,j,i) - sa_init(k) )
	1846	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
	1847	!
	1848	!-- Tendencies for the next Runge-Kutta step
	1849	IF ( runge_step == 1 ) THEN
	1850	tsa_m(k,j,i) = tend(k,j,i)
	1851	ELSEIF ( runge_step == 2 ) THEN
	1852	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tsa_m(k,j,i)
	1853	ENDIF
	1854	ENDDO
	1855	ENDDO
	1856	ENDDO
	1857
	1858	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
	1859
	1860	!
	1861	!-- Calculate density by the equation of state for seawater
	1862	CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
	1863	CALL eqn_state_seawater
	1864	CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
	1865
	1866	ENDIF
	1867
	1868	!
	1869	!-- If required, compute prognostic equation for total water content / scalar
	1870	IF ( humidity .OR. passive_scalar ) THEN
	1871
	1872	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
	1873
	1874	!
	1875	!-- Scalar/q-tendency terms with communication
	1876	sbt = tsc(2)
	1877	IF ( scalar_advec == 'bc-scheme' ) THEN
	1878
	1879	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1880	!
	1881	!-- Bott-Chlond scheme always uses Euler time step. Thus:
	1882	sbt = 1.0
	1883	ENDIF
	1884	tend = 0.0
	1885	CALL advec_s_bc( q, 'q' )
	1886
	1887	ENDIF
	1888
	1889	!
	1890	!-- Scalar/q-tendency terms with no communication
	1891	IF ( scalar_advec /= 'bc-scheme' ) THEN
	1892	tend = 0.0
	1893	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1894	IF ( ws_scheme_sca ) THEN
	1895	CALL advec_s_ws( q, 'q' )
	1896	ELSE
	1897	CALL advec_s_pw( q )
	1898	ENDIF
	1899	ELSE
	1900	CALL advec_s_up( q )
	1901	ENDIF
	1902	ENDIF
	1903
	1904	CALL diffusion_s( q, qsws, qswst, wall_qflux )
	1905
	1906	!
	1907	!-- If required compute decrease of total water content due to
	1908	!-- precipitation
	1909	IF ( precipitation ) THEN
	1910	CALL calc_precipitation
	1911	ENDIF
	1912
	1913	!
	1914	!-- Sink or source of scalar concentration due to canopy elements
	1915	IF ( plant_canopy ) CALL plant_canopy_model( 5 )
	1916
	1917	!
	1918	!-- If required compute influence of large-scale subsidence/ascent
	1919	IF ( large_scale_subsidence ) THEN
	1920	CALL subsidence( tend, q, q_init )
	1921	ENDIF
	1922
	1923	CALL user_actions( 'q-tendency' )
	1924
	1925	!
	1926	!-- Prognostic equation for total water content / scalar
	1927	DO i = nxl, nxr
	1928	DO j = nys, nyn
	1929	DO k = nzb_s_inner(j,i)+1, nzt
	1930	q_p(k,j,i) = q(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	1931	tsc(3) * tq_m(k,j,i) ) &
	1932	- tsc(5) * rdf_sc(k) * &
	1933	( q(k,j,i) - q_init(k) )
	1934	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
	1935	!
	1936	!-- Tendencies for the next Runge-Kutta step
	1937	IF ( runge_step == 1 ) THEN
	1938	tq_m(k,j,i) = tend(k,j,i)
	1939	ELSEIF ( runge_step == 2 ) THEN
	1940	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
	1941	ENDIF
	1942	ENDDO
	1943	ENDDO
	1944	ENDDO
	1945
	1946	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
	1947
	1948	ENDIF
	1949
	1950	!
	1951	!-- If required, compute prognostic equation for turbulent kinetic
	1952	!-- energy (TKE)
	1953	IF ( .NOT. constant_diffusion ) THEN
	1954
	1955	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
	1956
	1957	!
	1958	!-- TKE-tendency terms with communication
	1959	CALL production_e_init
	1960
	1961	sbt = tsc(2)
	1962	IF ( .NOT. use_upstream_for_tke ) THEN
	1963	IF ( scalar_advec == 'bc-scheme' ) THEN
	1964
	1965	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	1966	!
	1967	!-- Bott-Chlond scheme always uses Euler time step. Thus:
	1968	sbt = 1.0
	1969	ENDIF
	1970	tend = 0.0
	1971	CALL advec_s_bc( e, 'e' )
	1972
	1973	ENDIF
	1974	ENDIF
	1975
	1976	!
	1977	!-- TKE-tendency terms with no communication
	1978	IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN
	1979	IF ( use_upstream_for_tke ) THEN
	1980	tend = 0.0
	1981	CALL advec_s_up( e )
	1982	ELSE
	1983	IF ( timestep_scheme(1:5) == 'runge' ) THEN
	1984	IF ( ws_scheme_sca ) THEN
	1985	CALL advec_s_ws_acc( e, 'e' )
	1986	ELSE
	1987	tend = 0.0 ! to be removed later??
	1988	CALL advec_s_pw( e )
	1989	ENDIF
	1990	ELSE
	1991	tend = 0.0 ! to be removed later??
	1992	CALL advec_s_up( e )
	1993	ENDIF
	1994	ENDIF
	1995	ENDIF
	1996
	1997	IF ( .NOT. humidity ) THEN
	1998	IF ( ocean ) THEN
	1999	CALL diffusion_e( prho, prho_reference )
	2000	ELSE
	2001	CALL diffusion_e_acc( pt, pt_reference )
	2002	ENDIF
	2003	ELSE
	2004	CALL diffusion_e( vpt, pt_reference )
	2005	ENDIF
	2006
	2007	CALL production_e_acc
	2008
	2009	!
	2010	!-- Additional sink term for flows through plant canopies
	2011	IF ( plant_canopy ) CALL plant_canopy_model( 6 )
	2012	CALL user_actions( 'e-tendency' )
	2013
	2014	!
	2015	!-- Prognostic equation for TKE.
	2016	!-- Eliminate negative TKE values, which can occur due to numerical
	2017	!-- reasons in the course of the integration. In such cases the old TKE
	2018	!-- value is reduced by 90%.
	2019	!$acc kernels present( e, e_p, nzb_s_inner, tend, te_m )
	2020	!$acc loop
	2021	DO i = nxl, nxr
	2022	DO j = nys, nyn
	2023	!$acc loop vector( 32 )
	2024	DO k = 1, nzt
	2025	IF ( k > nzb_s_inner(j,i) ) THEN
	2026	e_p(k,j,i) = e(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
	2027	tsc(3) * te_m(k,j,i) )
	2028	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
	2029	!
	2030	!-- Tendencies for the next Runge-Kutta step
	2031	IF ( runge_step == 1 ) THEN
	2032	te_m(k,j,i) = tend(k,j,i)
	2033	ELSEIF ( runge_step == 2 ) THEN
	2034	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
	2035	ENDIF
	2036	ENDIF
	2037	ENDDO
	2038	ENDDO
	2039	ENDDO
	2040	!$acc end kernels
	2041
	2042	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
	2043	!$acc update host( e_p )
	2044
	2045	ENDIF
	2046
	2047
	2048	END SUBROUTINE prognostic_equations_acc
	2049
	2050
[736]	2051	END MODULE prognostic_equations_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |