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

Last change on this file since 2012 was 2012, checked in by kanani, 5 years ago

last commit documented

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