Home

Context Navigation

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

Last change on this file since 736 was 736, checked in by suehring, 13 years ago
Bugfix in ws-scheme concerning OpenMP paralellization.
Property svn:keywords set to `Id`
File size: 78.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |