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

Last change on this file since 1960 was 1960, checked in by suehring, 5 years ago

Separate balance equations for humidity and passive_scalar

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