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

Last change on this file since 1808 was 1758, checked in by maronga, 8 years ago

last commit documented

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