Home

Context Navigation

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

Last change on this file since 50 was 39, checked in by raasch, 18 years ago
comments prepared for 3.1c
Property svn:keywords set to `Id`
File size: 53.5 KB

Line
1	MODULE prognostic_equations_mod
2
3	!------------------------------------------------------------------------------!
4	! Actual revisions:
5	! -----------------
6	!
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: prognostic_equations.f90 39 2007-03-01 12:46:59Z raasch $
11	!
12	! 19 2007-02-23 04:53:48Z raasch
13	! Calculation of e, q, and pt extended for gridpoint nzt,
14	! handling of given temperature/humidity/scalar fluxes at top surface
15	!
16	! RCS Log replace by Id keyword, revision history cleaned up
17	!
18	! Revision 1.21 2006/08/04 15:01:07 raasch
19	! upstream scheme can be forced to be used for tke (use_upstream_for_tke)
20	! regardless of the timestep scheme used for the other quantities,
21	! new argument diss in call of diffusion_e
22	!
23	! Revision 1.1 2000/04/13 14:56:27 schroeter
24	! Initial revision
25	!
26	!
27	! Description:
28	! ------------
29	! Solving the prognostic equations.
30	!------------------------------------------------------------------------------!
31
32	USE arrays_3d
33	USE control_parameters
34	USE cpulog
35	USE grid_variables
36	USE indices
37	USE interfaces
38	USE pegrid
39	USE pointer_interfaces
40	USE statistics
41
42	USE advec_s_pw_mod
43	USE advec_s_up_mod
44	USE advec_u_pw_mod
45	USE advec_u_up_mod
46	USE advec_v_pw_mod
47	USE advec_v_up_mod
48	USE advec_w_pw_mod
49	USE advec_w_up_mod
50	USE buoyancy_mod
51	USE calc_precipitation_mod
52	USE calc_radiation_mod
53	USE coriolis_mod
54	USE diffusion_e_mod
55	USE diffusion_s_mod
56	USE diffusion_u_mod
57	USE diffusion_v_mod
58	USE diffusion_w_mod
59	USE impact_of_latent_heat_mod
60	USE production_e_mod
61	USE user_actions_mod
62
63
64	PRIVATE
65	PUBLIC prognostic_equations, prognostic_equations_fast, &
66	prognostic_equations_vec
67
68	INTERFACE prognostic_equations
69	MODULE PROCEDURE prognostic_equations
70	END INTERFACE prognostic_equations
71
72	INTERFACE prognostic_equations_fast
73	MODULE PROCEDURE prognostic_equations_fast
74	END INTERFACE prognostic_equations_fast
75
76	INTERFACE prognostic_equations_vec
77	MODULE PROCEDURE prognostic_equations_vec
78	END INTERFACE prognostic_equations_vec
79
80
81	CONTAINS
82
83
84	SUBROUTINE prognostic_equations
85
86	!------------------------------------------------------------------------------!
87	! Version with single loop optimization
88	!
89	! (Optimized over each single prognostic equation.)
90	!------------------------------------------------------------------------------!
91
92	IMPLICIT NONE
93
94	CHARACTER (LEN=9) :: time_to_string
95	INTEGER :: i, j, k
96	REAL :: sat, sbt
97
98	!
99	!-- Calculate those variables needed in the tendency terms which need
100	!-- global communication
101	CALL calc_mean_pt_profile( pt, 4 )
102	IF ( moisture ) CALL calc_mean_pt_profile( vpt, 44 )
103
104	!
105	!-- u-velocity component
106	CALL cpu_log( log_point(5), 'u-equation', 'start' )
107
108	!
109	!-- u-tendency terms with communication
110	IF ( momentum_advec == 'ups-scheme' ) THEN
111	tend = 0.0
112	CALL advec_u_ups
113	ENDIF
114
115	!
116	!-- u-tendency terms with no communication
117	DO i = nxl, nxr+uxrp
118	DO j = nys, nyn
119	!
120	!-- Tendency terms
121	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
122	tend(:,j,i) = 0.0
123	CALL advec_u_pw( i, j )
124	ELSE
125	IF ( momentum_advec /= 'ups-scheme' ) THEN
126	tend(:,j,i) = 0.0
127	CALL advec_u_up( i, j )
128	ENDIF
129	ENDIF
130	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
131	CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, u_m, &
132	usws_m, v_m, w_m, z0 )
133	ELSE
134	CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
135	v, w, z0 )
136	ENDIF
137	CALL coriolis( i, j, 1 )
138	IF ( sloping_surface ) CALL buoyancy( i, j, pt, 1, 4 )
139	CALL user_actions( i, j, 'u-tendency' )
140
141	!
142	!-- Prognostic equation for u-velocity component
143	DO k = nzb_u_inner(j,i)+1, nzt
144	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
145	dt_3d * ( &
146	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
147	- tsc(4) * ( p(k,j,i) - p(k,j,i-1) ) * ddx &
148	) - &
149	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
150	ENDDO
151
152	!
153	!-- Calculate tendencies for the next Runge-Kutta step
154	IF ( timestep_scheme(1:5) == 'runge' ) THEN
155	IF ( intermediate_timestep_count == 1 ) THEN
156	DO k = nzb_u_inner(j,i)+1, nzt
157	tu_m(k,j,i) = tend(k,j,i)
158	ENDDO
159	ELSEIF ( intermediate_timestep_count < &
160	intermediate_timestep_count_max ) THEN
161	DO k = nzb_u_inner(j,i)+1, nzt
162	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
163	ENDDO
164	ENDIF
165	ENDIF
166
167	ENDDO
168	ENDDO
169
170	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
171
172	!
173	!-- v-velocity component
174	CALL cpu_log( log_point(6), 'v-equation', 'start' )
175
176	!
177	!-- v-tendency terms with communication
178	IF ( momentum_advec == 'ups-scheme' ) THEN
179	tend = 0.0
180	CALL advec_v_ups
181	ENDIF
182
183	!
184	!-- v-tendency terms with no communication
185	DO i = nxl, nxr
186	DO j = nys, nyn+vynp
187	!
188	!-- Tendency terms
189	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
190	tend(:,j,i) = 0.0
191	CALL advec_v_pw( i, j )
192	ELSE
193	IF ( momentum_advec /= 'ups-scheme' ) THEN
194	tend(:,j,i) = 0.0
195	CALL advec_v_up( i, j )
196	ENDIF
197	ENDIF
198	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
199	CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, u_m, &
200	v_m, vsws_m, w_m, z0 )
201	ELSE
202	CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
203	vsws, w, z0 )
204	ENDIF
205	CALL coriolis( i, j, 2 )
206	CALL user_actions( i, j, 'v-tendency' )
207
208	!
209	!-- Prognostic equation for v-velocity component
210	DO k = nzb_v_inner(j,i)+1, nzt
211	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
212	dt_3d * ( &
213	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
214	- tsc(4) * ( p(k,j,i) - p(k,j-1,i) ) * ddy &
215	) - &
216	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
217	ENDDO
218
219	!
220	!-- Calculate tendencies for the next Runge-Kutta step
221	IF ( timestep_scheme(1:5) == 'runge' ) THEN
222	IF ( intermediate_timestep_count == 1 ) THEN
223	DO k = nzb_v_inner(j,i)+1, nzt
224	tv_m(k,j,i) = tend(k,j,i)
225	ENDDO
226	ELSEIF ( intermediate_timestep_count < &
227	intermediate_timestep_count_max ) THEN
228	DO k = nzb_v_inner(j,i)+1, nzt
229	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
230	ENDDO
231	ENDIF
232	ENDIF
233
234	ENDDO
235	ENDDO
236
237	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
238
239	!
240	!-- w-velocity component
241	CALL cpu_log( log_point(7), 'w-equation', 'start' )
242
243	!
244	!-- w-tendency terms with communication
245	IF ( momentum_advec == 'ups-scheme' ) THEN
246	tend = 0.0
247	CALL advec_w_ups
248	ENDIF
249
250	!
251	!-- w-tendency terms with no communication
252	DO i = nxl, nxr
253	DO j = nys, nyn
254	!
255	!-- Tendency terms
256	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
257	tend(:,j,i) = 0.0
258	CALL advec_w_pw( i, j )
259	ELSE
260	IF ( momentum_advec /= 'ups-scheme' ) THEN
261	tend(:,j,i) = 0.0
262	CALL advec_w_up( i, j )
263	ENDIF
264	ENDIF
265	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
266	CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, km_damp_y, &
267	tend, u_m, v_m, w_m, z0 )
268	ELSE
269	CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
270	tend, u, v, w, z0 )
271	ENDIF
272	CALL coriolis( i, j, 3 )
273	IF ( .NOT. moisture ) THEN
274	CALL buoyancy( i, j, pt, 3, 4 )
275	ELSE
276	CALL buoyancy( i, j, vpt, 3, 44 )
277	ENDIF
278	CALL user_actions( i, j, 'w-tendency' )
279
280	!
281	!-- Prognostic equation for w-velocity component
282	DO k = nzb_w_inner(j,i)+1, nzt-1
283	w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
284	dt_3d * ( &
285	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
286	- tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
287	) - &
288	tsc(5) * rdf(k) * w(k,j,i)
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_w_inner(j,i)+1, nzt-1
296	tw_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_w_inner(j,i)+1, nzt-1
301	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
302	ENDDO
303	ENDIF
304	ENDIF
305
306	ENDDO
307	ENDDO
308
309	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
310
311	!
312	!-- potential temperature
313	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
314
315	!
316	!-- pt-tendency terms with communication
317	IF ( scalar_advec == 'bc-scheme' ) THEN
318	!
319	!-- Bott-Chlond scheme always uses Euler time step. Thus:
320	sat = 1.0
321	sbt = 1.0
322	tend = 0.0
323	CALL advec_s_bc( pt, 'pt' )
324	ELSE
325	sat = tsc(1)
326	sbt = tsc(2)
327	IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN
328	tend = 0.0
329	CALL advec_s_ups( pt, 'pt' )
330	ENDIF
331	ENDIF
332
333	!
334	!-- pt-tendency terms with no communication
335	DO i = nxl, nxr
336	DO j = nys, nyn
337	!
338	!-- Tendency terms
339	IF ( scalar_advec == 'bc-scheme' ) THEN
340	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
341	ELSE
342	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
343	tend(:,j,i) = 0.0
344	CALL advec_s_pw( i, j, pt )
345	ELSE
346	IF ( scalar_advec /= 'ups-scheme' ) THEN
347	tend(:,j,i) = 0.0
348	CALL advec_s_up( i, j, pt )
349	ENDIF
350	ENDIF
351	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
352	THEN
353	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
354	tswst_m, tend )
355	ELSE
356	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
357	ENDIF
358	ENDIF
359
360	!
361	!-- If required compute heating/cooling due to long wave radiation
362	!-- processes
363	IF ( radiation ) THEN
364	CALL calc_radiation( i, j )
365	ENDIF
366
367	!
368	!-- If required compute impact of latent heat due to precipitation
369	IF ( precipitation ) THEN
370	CALL impact_of_latent_heat( i, j )
371	ENDIF
372	CALL user_actions( i, j, 'pt-tendency' )
373
374	!
375	!-- Prognostic equation for potential temperature
376	DO k = nzb_s_inner(j,i)+1, nzt
377	pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
378	dt_3d * ( &
379	sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
380	) - &
381	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
382	ENDDO
383
384	!
385	!-- Calculate tendencies for the next Runge-Kutta step
386	IF ( timestep_scheme(1:5) == 'runge' ) THEN
387	IF ( intermediate_timestep_count == 1 ) THEN
388	DO k = nzb_s_inner(j,i)+1, nzt
389	tpt_m(k,j,i) = tend(k,j,i)
390	ENDDO
391	ELSEIF ( intermediate_timestep_count < &
392	intermediate_timestep_count_max ) THEN
393	DO k = nzb_s_inner(j,i)+1, nzt
394	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
395	ENDDO
396	ENDIF
397	ENDIF
398
399	ENDDO
400	ENDDO
401
402	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
403
404	!
405	!-- If required, compute prognostic equation for total water content / scalar
406	IF ( moisture .OR. passive_scalar ) THEN
407
408	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
409
410	!
411	!-- Scalar/q-tendency terms with communication
412	IF ( scalar_advec == 'bc-scheme' ) THEN
413	!
414	!-- Bott-Chlond scheme always uses Euler time step. Thus:
415	sat = 1.0
416	sbt = 1.0
417	tend = 0.0
418	CALL advec_s_bc( q, 'q' )
419	ELSE
420	sat = tsc(1)
421	sbt = tsc(2)
422	IF ( tsc(2) /= 2.0 ) THEN
423	IF ( scalar_advec == 'ups-scheme' ) THEN
424	tend = 0.0
425	CALL advec_s_ups( q, 'q' )
426	ENDIF
427	ENDIF
428	ENDIF
429
430	!
431	!-- Scalar/q-tendency terms with no communication
432	DO i = nxl, nxr
433	DO j = nys, nyn
434	!
435	!-- Tendency-terms
436	IF ( scalar_advec == 'bc-scheme' ) THEN
437	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, tend )
438	ELSE
439	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
440	tend(:,j,i) = 0.0
441	CALL advec_s_pw( i, j, q )
442	ELSE
443	IF ( scalar_advec /= 'ups-scheme' ) THEN
444	tend(:,j,i) = 0.0
445	CALL advec_s_up( i, j, q )
446	ENDIF
447	ENDIF
448	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
449	THEN
450	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
451	qswst_m, tend )
452	ELSE
453	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
454	tend )
455	ENDIF
456	ENDIF
457
458	!
459	!-- If required compute decrease of total water content due to
460	!-- precipitation
461	IF ( precipitation ) THEN
462	CALL calc_precipitation( i, j )
463	ENDIF
464	CALL user_actions( i, j, 'q-tendency' )
465
466	!
467	!-- Prognostic equation for total water content / scalar
468	DO k = nzb_s_inner(j,i)+1, nzt
469	q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + &
470	dt_3d * ( &
471	sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
472	) - &
473	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
474	ENDDO
475
476	!
477	!-- Calculate tendencies for the next Runge-Kutta step
478	IF ( timestep_scheme(1:5) == 'runge' ) THEN
479	IF ( intermediate_timestep_count == 1 ) THEN
480	DO k = nzb_s_inner(j,i)+1, nzt
481	tq_m(k,j,i) = tend(k,j,i)
482	ENDDO
483	ELSEIF ( intermediate_timestep_count < &
484	intermediate_timestep_count_max ) THEN
485	DO k = nzb_s_inner(j,i)+1, nzt
486	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
487	ENDDO
488	ENDIF
489	ENDIF
490
491	ENDDO
492	ENDDO
493
494	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
495
496	ENDIF
497
498	!
499	!-- If required, compute prognostic equation for turbulent kinetic
500	!-- energy (TKE)
501	IF ( .NOT. constant_diffusion ) THEN
502
503	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
504
505	!
506	!-- TKE-tendency terms with communication
507	CALL production_e_init
508	IF ( .NOT. use_upstream_for_tke ) THEN
509	IF ( scalar_advec == 'bc-scheme' ) THEN
510	!
511	!-- Bott-Chlond scheme always uses Euler time step. Thus:
512	sat = 1.0
513	sbt = 1.0
514	tend = 0.0
515	CALL advec_s_bc( e, 'e' )
516	ELSE
517	sat = tsc(1)
518	sbt = tsc(2)
519	IF ( tsc(2) /= 2.0 ) THEN
520	IF ( scalar_advec == 'ups-scheme' ) THEN
521	tend = 0.0
522	CALL advec_s_ups( e, 'e' )
523	ENDIF
524	ENDIF
525	ENDIF
526	ENDIF
527
528	!
529	!-- TKE-tendency terms with no communication
530	DO i = nxl, nxr
531	DO j = nys, nyn
532	!
533	!-- Tendency-terms
534	IF ( scalar_advec == 'bc-scheme' .AND. &
535	.NOT. use_upstream_for_tke ) THEN
536	IF ( .NOT. moisture ) THEN
537	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
538	l_grid, pt, rif, tend, zu )
539	ELSE
540	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
541	l_grid, vpt, rif, tend, zu )
542	ENDIF
543	ELSE
544	IF ( use_upstream_for_tke ) THEN
545	tend(:,j,i) = 0.0
546	CALL advec_s_up( i, j, e )
547	ELSE
548	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
549	THEN
550	tend(:,j,i) = 0.0
551	CALL advec_s_pw( i, j, e )
552	ELSE
553	IF ( scalar_advec /= 'ups-scheme' ) THEN
554	tend(:,j,i) = 0.0
555	CALL advec_s_up( i, j, e )
556	ENDIF
557	ENDIF
558	ENDIF
559	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
560	THEN
561	IF ( .NOT. moisture ) THEN
562	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
563	km_m, l_grid, pt_m, rif_m, tend, zu )
564	ELSE
565	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
566	km_m, l_grid, vpt_m, rif_m, tend, zu )
567	ENDIF
568	ELSE
569	IF ( .NOT. moisture ) THEN
570	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
571	l_grid, pt, rif, tend, zu )
572	ELSE
573	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
574	l_grid, vpt, rif, tend, zu )
575	ENDIF
576	ENDIF
577	ENDIF
578	CALL production_e( i, j )
579	CALL user_actions( i, j, 'e-tendency' )
580
581	!
582	!-- Prognostic equation for TKE.
583	!-- Eliminate negative TKE values, which can occur due to numerical
584	!-- reasons in the course of the integration. In such cases the old TKE
585	!-- value is reduced by 90%.
586	DO k = nzb_s_inner(j,i)+1, nzt
587	e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + &
588	dt_3d * ( &
589	sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
590	)
591	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
592	ENDDO
593
594	!
595	!-- Calculate tendencies for the next Runge-Kutta step
596	IF ( timestep_scheme(1:5) == 'runge' ) THEN
597	IF ( intermediate_timestep_count == 1 ) THEN
598	DO k = nzb_s_inner(j,i)+1, nzt
599	te_m(k,j,i) = tend(k,j,i)
600	ENDDO
601	ELSEIF ( intermediate_timestep_count < &
602	intermediate_timestep_count_max ) THEN
603	DO k = nzb_s_inner(j,i)+1, nzt
604	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
605	ENDDO
606	ENDIF
607	ENDIF
608
609	ENDDO
610	ENDDO
611
612	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
613
614	ENDIF
615
616
617	END SUBROUTINE prognostic_equations
618
619
620	SUBROUTINE prognostic_equations_fast
621
622	!------------------------------------------------------------------------------!
623	! Version with one optimized loop over all equations. It is only allowed to
624	! be called for the standard Piascek-Williams advection scheme.
625	!
626	! The call of this subroutine is embedded in two DO loops over i and j, thus
627	! communication between CPUs is not allowed in this subroutine.
628	!
629	! (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
630	!------------------------------------------------------------------------------!
631
632	IMPLICIT NONE
633
634	CHARACTER (LEN=9) :: time_to_string
635	INTEGER :: i, j, k
636
637
638	!
639	!-- Time measurement can only be performed for the whole set of equations
640	CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
641
642
643	!
644	!-- Calculate those variables needed in the tendency terms which need
645	!-- global communication
646	CALL calc_mean_pt_profile( pt, 4 )
647	IF ( moisture ) CALL calc_mean_pt_profile( vpt, 44 )
648	IF ( .NOT. constant_diffusion ) CALL production_e_init
649
650
651	!
652	!-- Loop over all prognostic equations
653	!$OMP PARALLEL private (i,j,k)
654	!$OMP DO
655	DO i = nxl, nxr+uxrp ! Additional levels for non cyclic boundary
656	DO j = nys, nyn+vynp ! conditions are included
657	!
658	!-- Tendency terms for u-velocity component
659	IF ( j < nyn+1 ) THEN
660
661	tend(:,j,i) = 0.0
662	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
663	CALL advec_u_pw( i, j )
664	ELSE
665	CALL advec_u_up( i, j )
666	ENDIF
667	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
668	THEN
669	CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, &
670	u_m, usws_m, v_m, w_m, z0 )
671	ELSE
672	CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, &
673	usws, v, w, z0 )
674	ENDIF
675	CALL coriolis( i, j, 1 )
676	IF ( sloping_surface ) CALL buoyancy( i, j, pt, 1, 4 )
677	CALL user_actions( i, j, 'u-tendency' )
678
679	!
680	!-- Prognostic equation for u-velocity component
681	DO k = nzb_u_inner(j,i)+1, nzt
682	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
683	dt_3d * ( &
684	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
685	- tsc(4) * ( p(k,j,i) - p(k,j,i-1) ) * ddx &
686	) - &
687	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
688	ENDDO
689
690	!
691	!-- Calculate tendencies for the next Runge-Kutta step
692	IF ( timestep_scheme(1:5) == 'runge' ) THEN
693	IF ( intermediate_timestep_count == 1 ) THEN
694	DO k = nzb_u_inner(j,i)+1, nzt
695	tu_m(k,j,i) = tend(k,j,i)
696	ENDDO
697	ELSEIF ( intermediate_timestep_count < &
698	intermediate_timestep_count_max ) THEN
699	DO k = nzb_u_inner(j,i)+1, nzt
700	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
701	ENDDO
702	ENDIF
703	ENDIF
704
705	ENDIF
706
707	!
708	!-- Tendency terms for v-velocity component
709	IF ( i < nxr+1 ) THEN
710
711	tend(:,j,i) = 0.0
712	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
713	CALL advec_v_pw( i, j )
714	ELSE
715	CALL advec_v_up( i, j )
716	ENDIF
717	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
718	THEN
719	CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, &
720	u_m, v_m, vsws_m, w_m, z0 )
721	ELSE
722	CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
723	vsws, w, z0 )
724	ENDIF
725	CALL coriolis( i, j, 2 )
726	CALL user_actions( i, j, 'v-tendency' )
727
728	!
729	!-- Prognostic equation for v-velocity component
730	DO k = nzb_v_inner(j,i)+1, nzt
731	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
732	dt_3d * ( &
733	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
734	- tsc(4) * ( p(k,j,i) - p(k,j-1,i) ) * ddy &
735	) - &
736	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
737	ENDDO
738
739	!
740	!-- Calculate tendencies for the next Runge-Kutta step
741	IF ( timestep_scheme(1:5) == 'runge' ) THEN
742	IF ( intermediate_timestep_count == 1 ) THEN
743	DO k = nzb_v_inner(j,i)+1, nzt
744	tv_m(k,j,i) = tend(k,j,i)
745	ENDDO
746	ELSEIF ( intermediate_timestep_count < &
747	intermediate_timestep_count_max ) THEN
748	DO k = nzb_v_inner(j,i)+1, nzt
749	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
750	ENDDO
751	ENDIF
752	ENDIF
753
754	ENDIF
755
756	!
757	!-- Tendency terms for w-velocity component
758	IF ( i < nxr+1 .AND. j < nyn+1 ) THEN
759
760	tend(:,j,i) = 0.0
761	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
762	CALL advec_w_pw( i, j )
763	ELSE
764	CALL advec_w_up( i, j )
765	ENDIF
766	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
767	THEN
768	CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, &
769	km_damp_y, tend, u_m, v_m, w_m, z0 )
770	ELSE
771	CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
772	tend, u, v, w, z0 )
773	ENDIF
774	CALL coriolis( i, j, 3 )
775	IF ( .NOT. moisture ) THEN
776	CALL buoyancy( i, j, pt, 3, 4 )
777	ELSE
778	CALL buoyancy( i, j, vpt, 3, 44 )
779	ENDIF
780	CALL user_actions( i, j, 'w-tendency' )
781
782	!
783	!-- Prognostic equation for w-velocity component
784	DO k = nzb_w_inner(j,i)+1, nzt-1
785	w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
786	dt_3d * ( &
787	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
788	- tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
789	) - &
790	tsc(5) * rdf(k) * w(k,j,i)
791	ENDDO
792
793	!
794	!-- Calculate tendencies for the next Runge-Kutta step
795	IF ( timestep_scheme(1:5) == 'runge' ) THEN
796	IF ( intermediate_timestep_count == 1 ) THEN
797	DO k = nzb_w_inner(j,i)+1, nzt-1
798	tw_m(k,j,i) = tend(k,j,i)
799	ENDDO
800	ELSEIF ( intermediate_timestep_count < &
801	intermediate_timestep_count_max ) THEN
802	DO k = nzb_w_inner(j,i)+1, nzt-1
803	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
804	ENDDO
805	ENDIF
806	ENDIF
807
808	!
809	!-- Tendency terms for potential temperature
810	tend(:,j,i) = 0.0
811	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
812	CALL advec_s_pw( i, j, pt )
813	ELSE
814	CALL advec_s_up( i, j, pt )
815	ENDIF
816	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
817	THEN
818	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
819	tswst_m, tend )
820	ELSE
821	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
822	ENDIF
823
824	!
825	!-- If required compute heating/cooling due to long wave radiation
826	!-- processes
827	IF ( radiation ) THEN
828	CALL calc_radiation( i, j )
829	ENDIF
830
831	!
832	!-- If required compute impact of latent heat due to precipitation
833	IF ( precipitation ) THEN
834	CALL impact_of_latent_heat( i, j )
835	ENDIF
836	CALL user_actions( i, j, 'pt-tendency' )
837
838	!
839	!-- Prognostic equation for potential temperature
840	DO k = nzb_s_inner(j,i)+1, nzt
841	pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +&
842	dt_3d * ( &
843	tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
844	) - &
845	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
846	ENDDO
847
848	!
849	!-- Calculate tendencies for the next Runge-Kutta step
850	IF ( timestep_scheme(1:5) == 'runge' ) THEN
851	IF ( intermediate_timestep_count == 1 ) THEN
852	DO k = nzb_s_inner(j,i)+1, nzt
853	tpt_m(k,j,i) = tend(k,j,i)
854	ENDDO
855	ELSEIF ( intermediate_timestep_count < &
856	intermediate_timestep_count_max ) THEN
857	DO k = nzb_s_inner(j,i)+1, nzt
858	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
859	5.3125 * tpt_m(k,j,i)
860	ENDDO
861	ENDIF
862	ENDIF
863
864	!
865	!-- If required, compute prognostic equation for total water content /
866	!-- scalar
867	IF ( moisture .OR. passive_scalar ) THEN
868
869	!
870	!-- Tendency-terms for total water content / scalar
871	tend(:,j,i) = 0.0
872	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
873	THEN
874	CALL advec_s_pw( i, j, q )
875	ELSE
876	CALL advec_s_up( i, j, q )
877	ENDIF
878	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
879	THEN
880	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
881	qswst_m, tend )
882	ELSE
883	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
884	tend )
885	ENDIF
886
887	!
888	!-- If required compute decrease of total water content due to
889	!-- precipitation
890	IF ( precipitation ) THEN
891	CALL calc_precipitation( i, j )
892	ENDIF
893	CALL user_actions( i, j, 'q-tendency' )
894
895	!
896	!-- Prognostic equation for total water content / scalar
897	DO k = nzb_s_inner(j,i)+1, nzt
898	q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +&
899	dt_3d * ( &
900	tsc(2) * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
901	) - &
902	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
903	ENDDO
904
905	!
906	!-- Calculate tendencies for the next Runge-Kutta step
907	IF ( timestep_scheme(1:5) == 'runge' ) THEN
908	IF ( intermediate_timestep_count == 1 ) THEN
909	DO k = nzb_s_inner(j,i)+1, nzt
910	tq_m(k,j,i) = tend(k,j,i)
911	ENDDO
912	ELSEIF ( intermediate_timestep_count < &
913	intermediate_timestep_count_max ) THEN
914	DO k = nzb_s_inner(j,i)+1, nzt
915	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + &
916	5.3125 * tq_m(k,j,i)
917	ENDDO
918	ENDIF
919	ENDIF
920
921	ENDIF
922
923	!
924	!-- If required, compute prognostic equation for turbulent kinetic
925	!-- energy (TKE)
926	IF ( .NOT. constant_diffusion ) THEN
927
928	!
929	!-- Tendency-terms for TKE
930	tend(:,j,i) = 0.0
931	IF ( ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
932	.AND. .NOT. use_upstream_for_tke ) THEN
933	CALL advec_s_pw( i, j, e )
934	ELSE
935	CALL advec_s_up( i, j, e )
936	ENDIF
937	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
938	THEN
939	IF ( .NOT. moisture ) THEN
940	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
941	km_m, l_grid, pt_m, rif_m, tend, zu )
942	ELSE
943	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
944	km_m, l_grid, vpt_m, rif_m, tend, zu )
945	ENDIF
946	ELSE
947	IF ( .NOT. moisture ) THEN
948	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
949	l_grid, pt, rif, tend, zu )
950	ELSE
951	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
952	l_grid, vpt, rif, tend, zu )
953	ENDIF
954	ENDIF
955	CALL production_e( i, j )
956	CALL user_actions( i, j, 'e-tendency' )
957
958	!
959	!-- Prognostic equation for TKE.
960	!-- Eliminate negative TKE values, which can occur due to numerical
961	!-- reasons in the course of the integration. In such cases the old
962	!-- TKE value is reduced by 90%.
963	DO k = nzb_s_inner(j,i)+1, nzt
964	e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +&
965	dt_3d * ( &
966	tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
967	)
968	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
969	ENDDO
970
971	!
972	!-- Calculate tendencies for the next Runge-Kutta step
973	IF ( timestep_scheme(1:5) == 'runge' ) THEN
974	IF ( intermediate_timestep_count == 1 ) THEN
975	DO k = nzb_s_inner(j,i)+1, nzt
976	te_m(k,j,i) = tend(k,j,i)
977	ENDDO
978	ELSEIF ( intermediate_timestep_count < &
979	intermediate_timestep_count_max ) THEN
980	DO k = nzb_s_inner(j,i)+1, nzt
981	te_m(k,j,i) = -9.5625 * tend(k,j,i) + &
982	5.3125 * te_m(k,j,i)
983	ENDDO
984	ENDIF
985	ENDIF
986
987	ENDIF ! TKE equation
988
989	ENDIF ! Gridpoints excluding the non-cyclic wall
990
991	ENDDO
992	ENDDO
993	!$OMP END PARALLEL
994
995	CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
996
997
998	END SUBROUTINE prognostic_equations_fast
999
1000
1001	SUBROUTINE prognostic_equations_vec
1002
1003	!------------------------------------------------------------------------------!
1004	! Version for vector machines
1005	!------------------------------------------------------------------------------!
1006
1007	IMPLICIT NONE
1008
1009	CHARACTER (LEN=9) :: time_to_string
1010	INTEGER :: i, j, k
1011	REAL :: sat, sbt
1012
1013	!
1014	!-- Calculate those variables needed in the tendency terms which need
1015	!-- global communication
1016	CALL calc_mean_pt_profile( pt, 4 )
1017	IF ( moisture ) CALL calc_mean_pt_profile( vpt, 44 )
1018
1019	!
1020	!-- u-velocity component
1021	CALL cpu_log( log_point(5), 'u-equation', 'start' )
1022
1023	!
1024	!-- u-tendency terms with communication
1025	IF ( momentum_advec == 'ups-scheme' ) THEN
1026	tend = 0.0
1027	CALL advec_u_ups
1028	ENDIF
1029
1030	!
1031	!-- u-tendency terms with no communication
1032	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1033	tend = 0.0
1034	CALL advec_u_pw
1035	ELSE
1036	IF ( momentum_advec /= 'ups-scheme' ) THEN
1037	tend = 0.0
1038	CALL advec_u_up
1039	ENDIF
1040	ENDIF
1041	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1042	CALL diffusion_u( ddzu, ddzw, km_m, km_damp_y, tend, u_m, usws_m, v_m, &
1043	w_m, z0 )
1044	ELSE
1045	CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, w, z0 )
1046	ENDIF
1047	CALL coriolis( 1 )
1048	IF ( sloping_surface ) CALL buoyancy( pt, 1, 4 )
1049	CALL user_actions( 'u-tendency' )
1050
1051	!
1052	!-- Prognostic equation for u-velocity component
1053	DO i = nxl, nxr+uxrp
1054	DO j = nys, nyn
1055	DO k = nzb_u_inner(j,i)+1, nzt
1056	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
1057	dt_3d * ( &
1058	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
1059	- tsc(4) * ( p(k,j,i) - p(k,j,i-1) ) * ddx &
1060	) - &
1061	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
1062	ENDDO
1063	ENDDO
1064	ENDDO
1065
1066	!
1067	!-- Calculate tendencies for the next Runge-Kutta step
1068	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1069	IF ( intermediate_timestep_count == 1 ) THEN
1070	DO i = nxl, nxr+uxrp
1071	DO j = nys, nyn
1072	DO k = nzb_u_inner(j,i)+1, nzt
1073	tu_m(k,j,i) = tend(k,j,i)
1074	ENDDO
1075	ENDDO
1076	ENDDO
1077	ELSEIF ( intermediate_timestep_count < &
1078	intermediate_timestep_count_max ) THEN
1079	DO i = nxl, nxr+uxrp
1080	DO j = nys, nyn
1081	DO k = nzb_u_inner(j,i)+1, nzt
1082	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
1083	ENDDO
1084	ENDDO
1085	ENDDO
1086	ENDIF
1087	ENDIF
1088
1089	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1090
1091	!
1092	!-- v-velocity component
1093	CALL cpu_log( log_point(6), 'v-equation', 'start' )
1094
1095	!
1096	!-- v-tendency terms with communication
1097	IF ( momentum_advec == 'ups-scheme' ) THEN
1098	tend = 0.0
1099	CALL advec_v_ups
1100	ENDIF
1101
1102	!
1103	!-- v-tendency terms with no communication
1104	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1105	tend = 0.0
1106	CALL advec_v_pw
1107	ELSE
1108	IF ( momentum_advec /= 'ups-scheme' ) THEN
1109	tend = 0.0
1110	CALL advec_v_up
1111	ENDIF
1112	ENDIF
1113	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1114	CALL diffusion_v( ddzu, ddzw, km_m, km_damp_x, tend, u_m, v_m, vsws_m, &
1115	w_m, z0 )
1116	ELSE
1117	CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w, z0 )
1118	ENDIF
1119	CALL coriolis( 2 )
1120	CALL user_actions( 'v-tendency' )
1121
1122	!
1123	!-- Prognostic equation for v-velocity component
1124	DO i = nxl, nxr
1125	DO j = nys, nyn+vynp
1126	DO k = nzb_v_inner(j,i)+1, nzt
1127	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
1128	dt_3d * ( &
1129	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
1130	- tsc(4) * ( p(k,j,i) - p(k,j-1,i) ) * ddy &
1131	) - &
1132	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
1133	ENDDO
1134	ENDDO
1135	ENDDO
1136
1137	!
1138	!-- Calculate tendencies for the next Runge-Kutta step
1139	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1140	IF ( intermediate_timestep_count == 1 ) THEN
1141	DO i = nxl, nxr
1142	DO j = nys, nyn+vynp
1143	DO k = nzb_v_inner(j,i)+1, nzt
1144	tv_m(k,j,i) = tend(k,j,i)
1145	ENDDO
1146	ENDDO
1147	ENDDO
1148	ELSEIF ( intermediate_timestep_count < &
1149	intermediate_timestep_count_max ) THEN
1150	DO i = nxl, nxr
1151	DO j = nys, nyn+vynp
1152	DO k = nzb_v_inner(j,i)+1, nzt
1153	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
1154	ENDDO
1155	ENDDO
1156	ENDDO
1157	ENDIF
1158	ENDIF
1159
1160	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1161
1162	!
1163	!-- w-velocity component
1164	CALL cpu_log( log_point(7), 'w-equation', 'start' )
1165
1166	!
1167	!-- w-tendency terms with communication
1168	IF ( momentum_advec == 'ups-scheme' ) THEN
1169	tend = 0.0
1170	CALL advec_w_ups
1171	ENDIF
1172
1173	!
1174	!-- w-tendency terms with no communication
1175	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1176	tend = 0.0
1177	CALL advec_w_pw
1178	ELSE
1179	IF ( momentum_advec /= 'ups-scheme' ) THEN
1180	tend = 0.0
1181	CALL advec_w_up
1182	ENDIF
1183	ENDIF
1184	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1185	CALL diffusion_w( ddzu, ddzw, km_m, km_damp_x, km_damp_y, tend, u_m, &
1186	v_m, w_m, z0 )
1187	ELSE
1188	CALL diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, w, &
1189	z0 )
1190	ENDIF
1191	CALL coriolis( 3 )
1192	IF ( .NOT. moisture ) THEN
1193	CALL buoyancy( pt, 3, 4 )
1194	ELSE
1195	CALL buoyancy( vpt, 3, 44 )
1196	ENDIF
1197	CALL user_actions( 'w-tendency' )
1198
1199	!
1200	!-- Prognostic equation for w-velocity component
1201	DO i = nxl, nxr
1202	DO j = nys, nyn
1203	DO k = nzb_w_inner(j,i)+1, nzt-1
1204	w_p(k,j,i) = ( 1-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
1205	dt_3d * ( &
1206	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
1207	- tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
1208	) - &
1209	tsc(5) * rdf(k) * w(k,j,i)
1210	ENDDO
1211	ENDDO
1212	ENDDO
1213
1214	!
1215	!-- Calculate tendencies for the next Runge-Kutta step
1216	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1217	IF ( intermediate_timestep_count == 1 ) THEN
1218	DO i = nxl, nxr
1219	DO j = nys, nyn
1220	DO k = nzb_w_inner(j,i)+1, nzt-1
1221	tw_m(k,j,i) = tend(k,j,i)
1222	ENDDO
1223	ENDDO
1224	ENDDO
1225	ELSEIF ( intermediate_timestep_count < &
1226	intermediate_timestep_count_max ) THEN
1227	DO i = nxl, nxr
1228	DO j = nys, nyn
1229	DO k = nzb_w_inner(j,i)+1, nzt-1
1230	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
1231	ENDDO
1232	ENDDO
1233	ENDDO
1234	ENDIF
1235	ENDIF
1236
1237	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1238
1239	!
1240	!-- potential temperature
1241	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1242
1243	!
1244	!-- pt-tendency terms with communication
1245	IF ( scalar_advec == 'bc-scheme' ) THEN
1246	!
1247	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1248	sat = 1.0
1249	sbt = 1.0
1250	tend = 0.0
1251	CALL advec_s_bc( pt, 'pt' )
1252	ELSE
1253	sat = tsc(1)
1254	sbt = tsc(2)
1255	IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN
1256	tend = 0.0
1257	CALL advec_s_ups( pt, 'pt' )
1258	ENDIF
1259	ENDIF
1260
1261	!
1262	!-- pt-tendency terms with no communication
1263	IF ( scalar_advec == 'bc-scheme' ) THEN
1264	CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, tend )
1265	ELSE
1266	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1267	tend = 0.0
1268	CALL advec_s_pw( pt )
1269	ELSE
1270	IF ( scalar_advec /= 'ups-scheme' ) THEN
1271	tend = 0.0
1272	CALL advec_s_up( pt )
1273	ENDIF
1274	ENDIF
1275	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1276	CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tswst_m, tend )
1277	ELSE
1278	CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, tend )
1279	ENDIF
1280	ENDIF
1281
1282	!
1283	!-- If required compute heating/cooling due to long wave radiation
1284	!-- processes
1285	IF ( radiation ) THEN
1286	CALL calc_radiation
1287	ENDIF
1288
1289	!
1290	!-- If required compute impact of latent heat due to precipitation
1291	IF ( precipitation ) THEN
1292	CALL impact_of_latent_heat
1293	ENDIF
1294	CALL user_actions( 'pt-tendency' )
1295
1296	!
1297	!-- Prognostic equation for potential temperature
1298	DO i = nxl, nxr
1299	DO j = nys, nyn
1300	DO k = nzb_s_inner(j,i)+1, nzt
1301	pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
1302	dt_3d * ( &
1303	sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
1304	) - &
1305	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
1306	ENDDO
1307	ENDDO
1308	ENDDO
1309
1310	!
1311	!-- Calculate tendencies for the next Runge-Kutta step
1312	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1313	IF ( intermediate_timestep_count == 1 ) THEN
1314	DO i = nxl, nxr
1315	DO j = nys, nyn
1316	DO k = nzb_s_inner(j,i)+1, nzt
1317	tpt_m(k,j,i) = tend(k,j,i)
1318	ENDDO
1319	ENDDO
1320	ENDDO
1321	ELSEIF ( intermediate_timestep_count < &
1322	intermediate_timestep_count_max ) THEN
1323	DO i = nxl, nxr
1324	DO j = nys, nyn
1325	DO k = nzb_s_inner(j,i)+1, nzt
1326	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
1327	ENDDO
1328	ENDDO
1329	ENDDO
1330	ENDIF
1331	ENDIF
1332
1333	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1334
1335	!
1336	!-- If required, compute prognostic equation for total water content / scalar
1337	IF ( moisture .OR. passive_scalar ) THEN
1338
1339	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
1340
1341	!
1342	!-- Scalar/q-tendency terms with communication
1343	IF ( scalar_advec == 'bc-scheme' ) THEN
1344	!
1345	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1346	sat = 1.0
1347	sbt = 1.0
1348	tend = 0.0
1349	CALL advec_s_bc( q, 'q' )
1350	ELSE
1351	sat = tsc(1)
1352	sbt = tsc(2)
1353	IF ( tsc(2) /= 2.0 ) THEN
1354	IF ( scalar_advec == 'ups-scheme' ) THEN
1355	tend = 0.0
1356	CALL advec_s_ups( q, 'q' )
1357	ENDIF
1358	ENDIF
1359	ENDIF
1360
1361	!
1362	!-- Scalar/q-tendency terms with no communication
1363	IF ( scalar_advec == 'bc-scheme' ) THEN
1364	CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, tend )
1365	ELSE
1366	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1367	tend = 0.0
1368	CALL advec_s_pw( q )
1369	ELSE
1370	IF ( scalar_advec /= 'ups-scheme' ) THEN
1371	tend = 0.0
1372	CALL advec_s_up( q )
1373	ENDIF
1374	ENDIF
1375	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1376	CALL diffusion_s( ddzu, ddzw, kh_m, q_m, qsws_m, qswst_m, tend )
1377	ELSE
1378	CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, tend )
1379	ENDIF
1380	ENDIF
1381
1382	!
1383	!-- If required compute decrease of total water content due to
1384	!-- precipitation
1385	IF ( precipitation ) THEN
1386	CALL calc_precipitation
1387	ENDIF
1388	CALL user_actions( 'q-tendency' )
1389
1390	!
1391	!-- Prognostic equation for total water content / scalar
1392	DO i = nxl, nxr
1393	DO j = nys, nyn
1394	DO k = nzb_s_inner(j,i)+1, nzt
1395	q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + &
1396	dt_3d * ( &
1397	sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
1398	) - &
1399	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(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 = nys, nyn
1410	DO k = nzb_s_inner(j,i)+1, nzt
1411	tq_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 = nys, nyn
1419	DO k = nzb_s_inner(j,i)+1, nzt
1420	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
1421	ENDDO
1422	ENDDO
1423	ENDDO
1424	ENDIF
1425	ENDIF
1426
1427	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
1428
1429	ENDIF
1430
1431	!
1432	!-- If required, compute prognostic equation for turbulent kinetic
1433	!-- energy (TKE)
1434	IF ( .NOT. constant_diffusion ) THEN
1435
1436	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
1437
1438	!
1439	!-- TKE-tendency terms with communication
1440	CALL production_e_init
1441	IF ( .NOT. use_upstream_for_tke ) THEN
1442	IF ( scalar_advec == 'bc-scheme' ) THEN
1443	!
1444	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1445	sat = 1.0
1446	sbt = 1.0
1447	tend = 0.0
1448	CALL advec_s_bc( e, 'e' )
1449	ELSE
1450	sat = tsc(1)
1451	sbt = tsc(2)
1452	IF ( tsc(2) /= 2.0 ) THEN
1453	IF ( scalar_advec == 'ups-scheme' ) THEN
1454	tend = 0.0
1455	CALL advec_s_ups( e, 'e' )
1456	ENDIF
1457	ENDIF
1458	ENDIF
1459	ENDIF
1460
1461	!
1462	!-- TKE-tendency terms with no communication
1463	IF ( scalar_advec == 'bc-scheme' .AND. .NOT. use_upstream_for_tke ) &
1464	THEN
1465	IF ( .NOT. moisture ) THEN
1466	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
1467	rif, tend, zu )
1468	ELSE
1469	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
1470	rif, tend, zu )
1471	ENDIF
1472	ELSE
1473	IF ( use_upstream_for_tke ) THEN
1474	tend = 0.0
1475	CALL advec_s_up( e )
1476	ELSE
1477	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1478	tend = 0.0
1479	CALL advec_s_pw( e )
1480	ELSE
1481	IF ( scalar_advec /= 'ups-scheme' ) THEN
1482	tend = 0.0
1483	CALL advec_s_up( e )
1484	ENDIF
1485	ENDIF
1486	ENDIF
1487	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1488	IF ( .NOT. moisture ) THEN
1489	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
1490	pt_m, rif_m, tend, zu )
1491	ELSE
1492	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
1493	vpt_m, rif_m, tend, zu )
1494	ENDIF
1495	ELSE
1496	IF ( .NOT. moisture ) THEN
1497	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
1498	rif, tend, zu )
1499	ELSE
1500	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
1501	rif, tend, zu )
1502	ENDIF
1503	ENDIF
1504	ENDIF
1505	CALL production_e
1506	CALL user_actions( 'e-tendency' )
1507
1508	!
1509	!-- Prognostic equation for TKE.
1510	!-- Eliminate negative TKE values, which can occur due to numerical
1511	!-- reasons in the course of the integration. In such cases the old TKE
1512	!-- value is reduced by 90%.
1513	DO i = nxl, nxr
1514	DO j = nys, nyn
1515	DO k = nzb_s_inner(j,i)+1, nzt
1516	e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + &
1517	dt_3d * ( &
1518	sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
1519	)
1520	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
1521	ENDDO
1522	ENDDO
1523	ENDDO
1524
1525	!
1526	!-- Calculate tendencies for the next Runge-Kutta step
1527	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1528	IF ( intermediate_timestep_count == 1 ) THEN
1529	DO i = nxl, nxr
1530	DO j = nys, nyn
1531	DO k = nzb_s_inner(j,i)+1, nzt
1532	te_m(k,j,i) = tend(k,j,i)
1533	ENDDO
1534	ENDDO
1535	ENDDO
1536	ELSEIF ( intermediate_timestep_count < &
1537	intermediate_timestep_count_max ) THEN
1538	DO i = nxl, nxr
1539	DO j = nys, nyn
1540	DO k = nzb_s_inner(j,i)+1, nzt
1541	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
1542	ENDDO
1543	ENDDO
1544	ENDDO
1545	ENDIF
1546	ENDIF
1547
1548	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
1549
1550	ENDIF
1551
1552
1553	END SUBROUTINE prognostic_equations_vec
1554
1555
1556	END MODULE prognostic_equations_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |