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

Last change on this file since 1616 was 1586, checked in by maronga, 10 years ago

last commit documented

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