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

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