Home

Context Navigation

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

Last change on this file since 410 was 407, checked in by maronga, 14 years ago
humidity for non-flat topography implemented, re-adjustments for lcxt4
Property svn:keywords set to `Id`
File size: 70.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |