Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |