Home

Context Navigation

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

Last change on this file since 769 was 737, checked in by suehring, 13 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 78.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |