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

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

changes related to steering and formating of urban surface model

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