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

Last change on this file since 1374 was 1374, checked in by raasch, 7 years ago

bugfixes:
missing variables added to ONLY lists in USE statements (advec_s_bc, advec_s_pw, advec_s_up, advec_ws, buoyancy, diffusion_e, diffusion_s, fft_xy, flow_statistics, palm, prognostic_equations)
missing module kinds added (cuda_fft_interfaces)
dpk renamed dp (fft_xy)
missing dependency for check_open added (Makefile)
variables removed from acc-present-list (diffusion_e, diffusion_w, diffusivities, production_e, wall_fluxes)
syntax errors removed from openacc-branch (flow_statistics)
USE-statement for nopointer-case added (swap_timelevel)

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