Home

Context Navigation

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

Last change on this file since 1007 was 1002, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 62.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |