Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |