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

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