Home

Context Navigation

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

Last change on this file since 720 was 710, checked in by raasch, 14 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 75.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |