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

Last change on this file since 1412 was 1410, checked in by suehring, 10 years ago

last commit documented

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