Home

Context Navigation

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

Last change on this file since 940 was 940, checked in by raasch, 12 years ago
temperature equation can be switched off; bugfix of tridia_1dd for current Intel (12.1) compilers
Property svn:keywords set to `Id`
File size: 80.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |