Home

Context Navigation

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

Last change on this file since 785 was 785, checked in by raasch, 12 years ago
scalar quantities can be excluded from Rayleigh damping; bugfix for long lines in configuration file with more than 300 characters
Property svn:keywords set to `Id`
File size: 78.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |