Home

Context Navigation

prognostic_equations.f90 @ 404

Last change on this file since 404 was 392, checked in by raasch, 15 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:

lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:

Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list

Property svn:keywords set to Id

File size: 70.4 KB

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

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

Context Navigation

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

Download in other formats: