Home

Context Navigation

prognostic_equations.f90 @ 667

Last change on this file since 667 was 667, checked in by suehring, 13 years ago

summary:

Gryschka:

Coupling with different resolution and different numbers of PEs in ocean and atmosphere is available
Exchange of u and v from ocean surface to atmosphere surface
Mirror boundary condition for u and v at the bottom are replaced by dirichlet boundary conditions
Inflow turbulence is now defined by flucuations around spanwise mean
Bugfixes for cyclic_fill and constant_volume_flow

Suehring:

New advection added ( Wicker and Skamarock 5th order ), therefore:
- New module advec_ws.f90
- Modified exchange of ghost boundaries.
- Modified evaluation of turbulent fluxes
- New index bounds nxlg, nxrg, nysg, nyng

advec_ws.f90

Advection scheme for scalars and momentum using the flux formulation of
Wicker and Skamarock 5th order.
Additionally the module contains of a routine using for initialisation and
steering of the statical evaluation. The computation of turbulent fluxes takes
place inside the advection routines.
In case of vector architectures Dirichlet and Radiation boundary conditions are
outstanding and not available. Furthermore simulations within topography are
not possible so far. A further routine local_diss_ij is available and is used
if a control of dissipative fluxes is desired.

check_parameters.f90

Exchange of parameters between ocean and atmosphere via PE0
Check for illegal combination of ws-scheme and timestep scheme.
Check for topography and ws-scheme.
Check for not cyclic boundary conditions in combination with ws-scheme and
loop_optimization = 'vector'.
Check for call_psolver_at_all_substeps and ws-scheme for momentum_advec.

Different processor/grid topology in atmosphere and ocean is now allowed!
Bugfixes in checking for conserve_volume_flow_mode.

exchange_horiz.f90

Dynamic exchange of ghost points with nbgp_local to ensure that no useless
ghost points exchanged in case of multigrid. type_yz(0) and type_xz(0) used for
normal grid, the remaining types used for the several grid levels.
Exchange is done via MPI-Vectors with a dynamic value of ghost points which
depend on the advection scheme. Exchange of left and right PEs is 10% faster
with MPI-Vectors than without.

flow_statistics.f90

When advection is computed with ws-scheme, turbulent fluxes are already
computed in the respective advection routines and buffered in arrays
sums_xxxx_ws_l(). This is due to a consistent treatment of statistics
with the numerics and to avoid unphysical kinks near the surface. So some if-
requests has to be done to dicern between fluxes from ws-scheme other advection
schemes. Furthermore the computation of z_i is only done if the heat flux
exceeds a minimum value. This affects only simulations of a neutral boundary
layer and is due to reasons of computations in the advection scheme.

inflow_turbulence.f90

Using nbgp recycling planes for a better resolution of the turbulent flow near
the inflow.

init_grid.f90

Definition of new array bounds nxlg, nxrg, nysg, nyng on each PE.
Furthermore the allocation of arrays and steering of loops is done with these
parameters. Call of exchange_horiz are modified.
In case of dirichlet bounday condition at the bottom zu(0)=0.0
dzu_mg has to be set explicitly for a equally spaced grid near bottom.
ddzu_pres added to use a equally spaced grid near bottom.

init_pegrid.f90

Moved determination of target_id's from init_coupling
Determination of parameters needed for coupling (coupling_topology, ngp_a, ngp_o)
with different grid/processor-topology in ocean and atmosphere

Adaption of ngp_xy, ngp_y to a dynamic number of ghost points.
The maximum_grid_level changed from 1 to 0. 0 is the normal grid, 1 to
maximum_grid_level the grids for multigrid, in which 0 and 1 are normal grids.
This distinction is due to reasons of data exchange and performance for the
normal grid and grids in poismg.
The definition of MPI-Vectors adapted to a dynamic numer of ghost points.
New MPI-Vectors for data exchange between left and right boundaries added.
This is due to reasons of performance (10% faster).

ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
TEST OUTPUT (TO BE REMOVED) logging mpi2 ierr values

parin.f90

Steering parameter dissipation_control added in inipar.

Makefile

Module advec_ws added.

Modules

Removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc

For coupling with different resolution in ocean and atmophere:
+nx_a, +nx_o, ny_a, +ny_o, ngp_a, ngp_o, +total_2d_o, +total_2d_a,
+coupling_topology

Buffer arrays for the left sided advective fluxes added in arrays_3d.
+flux_s_u, +flux_s_v, +flux_s_w, +diss_s_u, +diss_s_v, +diss_s_w,
+flux_s_pt, +diss_s_pt, +flux_s_e, +diss_s_e, +flux_s_q, +diss_s_q,
+flux_s_sa, +diss_s_sa
3d arrays for dissipation control added. (only necessary for vector arch.)
+var_x, +var_y, +var_z, +gamma_x, +gamma_y, +gamma_z
Default of momentum_advec and scalar_advec changed to 'ws-scheme' .
+exchange_mg added in control_parameters to steer the data exchange.
Parameters +nbgp, +nxlg, +nxrg, +nysg, +nyng added in indices.
flag array +boundary_flags added in indices to steer the degradation of order
of the advective fluxes when non-cyclic boundaries are used.
MPI-datatypes +type_y, +type_y_int and +type_yz for data_exchange added in
pegrid.
+sums_wsus_ws_l, +sums_wsvs_ws_l, +sums_us2_ws_l, +sums_vs2_ws_l,
+sums_ws2_ws_l, +sums_wspts_ws_l, +sums_wssas_ws_l, +sums_wsqs_ws_l
and +weight_substep added in statistics to steer the statistical evaluation
of turbulent fluxes in the advection routines.
LOGICALS +ws_scheme_sca and +ws_scheme_mom added to get a better performance
in prognostic_equations.
LOGICAL +dissipation_control control added to steer numerical dissipation
in ws-scheme.

Changed length of string run_description_header

pres.f90

New allocation of tend when ws-scheme and multigrid is used. This is due to
reasons of perforance of the data_exchange. The same is done with p after
poismg is called.
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng when no
multigrid is used. Calls of exchange_horiz are modified.

bugfix: After pressure correction no volume flow correction in case of
non-cyclic boundary conditions
(has to be done only before pressure correction)

Call of SOR routine is referenced with ddzu_pres.

prognostic_equations.f90

Calls of the advection routines with WS5 added.
Calls of ws_statistics added to set the statistical arrays to zero after each
time step.

advec_particles.f90

Declaration of de_dx, de_dy, de_dz adapted to additional ghost points.
Furthermore the calls of exchange_horiz were modified.

asselin_filter.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

average_3d_data.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

boundary_conds.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
Removed mirror boundary conditions for u and v at the bottom in case of
ibc_uv_b == 0. Instead, dirichelt boundary conditions (u=v=0) are set
in init_3d_model

calc_liquid_water_content.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

calc_spectra.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for
allocation of tend.

check_open.f90

Output of total array size was adapted to nbgp.

data_output_2d.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays local_2d and total_2d.
Calls of exchange_horiz are modified.

data_output_2d.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Skip-value skip_do_avs changed to a dynamic adaption of ghost points.

data_output_mask.f90

Calls of exchange_horiz are modified.

diffusion_e.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_s.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_u.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_v.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_w.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

exchange_horiz_2d.f90

Dynamic exchange of ghost points with nbgp, which depends on the advection
scheme. Exchange between left and right PEs is now done with MPI-vectors.

global_min_max.f90

Adapting of the index arrays, because MINLOC assumes lowerbound
at 1 and not at nbgp.

init_3d_model.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Call ws_init to initialize arrays needed for statistical evaluation and
optimization when ws-scheme is used.
Initial volume flow is now calculated by using the variable hom_sum.
Therefore the correction of initial volume flow for non-flat topography
removed (removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc)
Changed surface boundary conditions for u and v in case of ibc_uv_b == 0 from
mirror bc to dirichlet boundary conditions (u=v=0), so that k=nzb is
representative for the height z0

Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner)

init_coupling.f90

determination of target_id's moved to init_pegrid

init_pt_anomaly.f90

Call of exchange_horiz are modified.

init_rankine.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

init_slope.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

header.f90

Output of advection scheme.

poismg.f90

Calls of exchange_horiz are modified.

prandtl_fluxes.f90

Changed surface boundary conditions for u and v from mirror bc to dirichelt bc,
therefore u(uzb,:,:) and v(nzb,:,:) is now representative for the height z0
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

production_e.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

read_3d_binary.f90

+/- 1 replaced with +/- nbgp when swapping and allocating variables.

sor.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Call of exchange_horiz are modified.
bug removed in declaration of ddzw(), nz replaced by nzt+1

subsidence.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

sum_up_3d_data.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

surface_coupler.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in
MPI_SEND() and MPI_RECV.
additional case for nonequivalent processor and grid topopolgy in ocean and
atmosphere added (coupling_topology = 1)

Added exchange of u and v from Ocean to Atmosphere

time_integration.f90

Calls of exchange_horiz are modified.
Adaption to slooping surface.

timestep.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_3d_data_averaging.f90, user_data_output_2d.f90, user_data_output_3d.f90,
user_actions.f90, user_init.f90, user_init_plant_canopy.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_read_restart_data.f90

Allocation with nbgp.

wall_fluxes.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

write_compressed.f90

Array bounds and nx, ny adapted with nbgp.

sor.f90

bug removed in declaration of ddzw(), nz replaced by nzt+1

Property svn:keywords set to Id

File size: 75.8 KB

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

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

Context Navigation

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

Download in other formats: