Home

Context Navigation

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

Last change on this file since 957 was 941, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 80.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |