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

Last change on this file since 1822 was 1822, checked in by hoffmann, 8 years ago

changes in LPM and bulk cloud microphysics

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