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

Last change on this file since 1518 was 1518, checked in by hoffmann, 7 years ago

last commit documented

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