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

Last change on this file since 2186 was 2156, checked in by hoffmann, 8 years ago

last commit documented

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