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

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

last commit documented

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