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

Last change on this file since 1823 was 1823, checked in by hoffmann, 5 years ago

last commit documented

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