Home

Context Navigation

prognostic_equations.f90 @ 1111

Last change on this file since 1111 was 1111, checked in by raasch, 11 years ago

New:
---

GPU porting of pres, swap_timelevel. Adjustments of openACC directives.
Further porting of poisfft, which now runs completely on GPU without any
host/device data transfer for serial an parallel runs (but parallel runs
require data transfer before and after the MPI transpositions).
GPU-porting of tridiagonal solver:
tridiagonal routines split into extermal subroutines (instead using CONTAINS),
no distinction between parallel/non-parallel in poisfft and tridia any more,
tridia routines moved to end of file because of probable bug in PGI compiler
(otherwise "invalid device function" is indicated during runtime).
(cuda_fft_interfaces, fft_xy, flow_statistics, init_3d_model, palm, poisfft, pres, prognostic_equations, swap_timelevel, time_integration, transpose)
output of accelerator board information. (header)

optimization of tridia routines: constant elements and coefficients of tri are
stored in seperate arrays ddzuw and tric, last dimension of tri reduced from 5 to 2,
(init_grid, init_3d_model, modules, palm, poisfft)

poisfft_init is now called internally from poisfft,
(Makefile, Makefile_check, init_pegrid, poisfft, poisfft_hybrid)

CPU-time per grid point and timestep is output to CPU_MEASURES file
(cpu_statistics, modules, time_integration)

Changed:

resorting from/to array work changed, work now has 4 dimensions instead of 1 (transpose)
array diss allocated only if required (init_3d_model)

pressure boundary condition "Neumann+inhomo" removed from the code
(check_parameters, header, poisfft, poisfft_hybrid, pres)

Errors:

bugfix: dependency added for cuda_fft_interfaces (Makefile)
bugfix: CUDA fft plans adjusted for domain decomposition (before they always
used total domain) (fft_xy)

Property svn:keywords set to Id

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

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

Context Navigation

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

Download in other formats: