Home

Context Navigation

prognostic_equations.f90 @ 80

Last change on this file since 80 was 77, checked in by raasch, 18 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:

General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:

Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

Property svn:keywords set to Id

File size: 54.6 KB

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

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

Context Navigation

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

Download in other formats: