Home

Context Navigation

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

Last change on this file since 709 was 709, checked in by raasch, 13 years ago
formatting adjustments
Property svn:keywords set to `Id`
File size: 75.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |