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

Last change on this file since 1874 was 1874, checked in by maronga, 5 years ago

last commit documented

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