Home

Context Navigation

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

Last change on this file since 978 was 978, checked in by fricke, 12 years ago
merge fricke branch back into trunk
Property svn:keywords set to `Id`
File size: 80.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |