Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |