Home

Context Navigation

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

Last change on this file since 992 was 979, checked in by fricke, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 80.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |