Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |