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

Last change on this file since 1933 was 1917, checked in by witha, 8 years ago

Corrected documented changes in the file headers

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