source: palm/trunk/SOURCE/prognostic_equations_mod.f90 @ 1850

Last change on this file since 1850 was 1850, checked in by maronga, 6 years ago

added _mod string to several filenames to meet the naming convection for modules

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