Home

Context Navigation

prognostic_equations.f90 @ 141

Last change on this file since 141 was 139, checked in by raasch, 17 years ago

New:
---

Plant canopy model of Watanabe (2004,BLM 112,307-341) added.
It can be switched on by the inipar parameter plant_canopy.
The inipar parameter canopy_mode can be used to prescribe a
plant canopy type. The default case is a homogeneous plant
canopy. Heterogeneous distributions of the leaf area
density and the canopy drag coefficient can be defined in the
new routine user_init_plant_canopy (user_interface).
The inipar parameters lad_surface, lad_vertical_gradient and
lad_vertical_gradient_level can be used in order to
prescribe the vertical profile of leaf area density. The
inipar parameter drag_coefficient determines the canopy
drag coefficient.
Finally, the inipar parameter pch_index determines the
index of the upper boundary of the plant canopy.

Allow new case bc_uv_t = 'dirichlet_0' for channel flow.

For unknown variables (CASE DEFAULT) call new subroutine user_data_output_dvrp

Pressure boundary conditions for vertical walls added to the multigrid solver.
They are applied using new wall flag arrays (wall_flags_..) which are defined
for each grid level. New argument gls added to routine user_init_grid
(user_interface).

Frequence of sorting particles can be controlled with new particles_par
parameter dt_sort_particles. Sorting is moved from the SGS timestep loop in
advec_particles after the end of this loop.

advec_particles, check_parameters, data_output_dvrp, header, init_3d_model, init_grid, init_particles, init_pegrid, modules, package_parin, parin, plant_canopy_model, read_var_list, read_3d_binary, user_interface, write_var_list, write_3d_binary

Changed:

Redefine initial nzb_local as the actual total size of topography (later the
extent of topography in nzb_local is reduced by 1dx at the E topography walls
and by 1dy at the N topography walls to form the basis for nzb_s_inner);
for consistency redefine 'single_building' case.

Vertical profiles now based on nzb_s_inner; they are divided by
ngp_2dh_s_inner (scalars, procucts of scalars) and ngp_2dh (staggered velocity
components and their products, procucts of scalars and velocity components),
respectively.

Allow two instead of one digit to specify isosurface and slicer variables.

Status of 3D-volume NetCDF data file only depends on switch netcdf_64bit_3d (check_open)

prognostic_equations include the respective wall_*flux in the parameter list of
calls of diffusion_s. Same as before, only the values of wall_heatflux(0:4)
can be assigned. At present, wall_humidityflux, wall_qflux, wall_salinityflux,
and wall_scalarflux are kept zero. diffusion_s uses the respective wall_*flux
instead of wall_heatflux. This update serves two purposes:

it avoids errors in calculations with humidity/scalar/salinity and prescribed

non-zero wall_heatflux,

it prepares PALM for a possible assignment of wall fluxes of

humidity/scalar/salinity in a future release.

buoyancy, check_open, data_output_dvrp, diffusion_s, diffusivities, flow_statistics, header, init_3d_model, init_dvrp, init_grid, modules, prognostic_equations

Errors:

Bugfix: summation of sums_l_l in diffusivities.

Several bugfixes in the ocean part: Initial density rho is calculated
(init_ocean). Error in initializing u_init and v_init removed
(check_parameters). Calculation of density flux now starts from
nzb+1 (production_e).

Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
numbers along y, small bugfixes in the SGS part (advec_particles)

Bugfix: model_string needed a default value (combine_plot_fields)

Bugfix: wavenumber calculation for even nx in routines maketri (poisfft)

Bugfix: assignment of fluxes at walls

Bugfix: absolute value of f must be used when calculating the Blackadar mixing length (init_1d_model)

advec_particles, check_parameters, combine_plot_fields, diffusion_s, diffusivities, init_ocean, init_1d_model, poisfft, production_e

Property svn:keywords set to Id

File size: 67.7 KB

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

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

Context Navigation

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

Download in other formats: