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

Last change on this file since 1496 was 1496, checked in by maronga, 9 years ago

added beta version of a land surface model and a simple radiation model for clear sky conditions

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