Home

Context Navigation

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

Last change on this file since 1016 was 1015, checked in by raasch, 12 years ago
Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added. Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed. mask array is set zero for ghost boundaries
Property svn:keywords set to `Id`
File size: 79.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |