Home

Context Navigation

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

Last change on this file since 420 was 411, checked in by heinze, 15 years ago
Large scale vertical motion (subsidence/ascent) can be applied to the prognostic equation for the potential temperature
Property svn:keywords set to `Id`
File size: 71.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |