Home

Context Navigation

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

Last change on this file since 809 was 786, checked in by raasch, 13 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 78.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |