Home

Context Navigation

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

Last change on this file since 2153 was 2119, checked in by raasch, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 66.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |