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

Last change on this file since 1685 was 1683, checked in by knoop, 9 years ago

last commit documented

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