Home

Context Navigation

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

Last change on this file since 203 was 198, checked in by raasch, 16 years ago
file headers updated for the next release 3.5
Property svn:keywords set to `Id`
File size: 68.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |