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

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

last commit documented

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