Home

Context Navigation

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

Last change on this file since 1120 was 1116, checked in by hoffmann, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 66.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |