Home

Context Navigation

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

Last change on this file since 673 was 673, checked in by suehring, 14 years ago
Right computation of the pressure using Runge-Kutta weighting coefficients. Consideration of the pressure gradient during the time integration removed. Removed bugfix concerning velocity variances.
Property svn:keywords set to `Id`
File size: 75.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |