Home

Context Navigation

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

Last change on this file since 1182 was 1182, checked in by raasch, 11 years ago
last commit documented, rc-file for example run updated
Property svn:keywords set to `Id`
File size: 65.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |