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

Last change on this file since 2192 was 2192, checked in by raasch, 5 years ago

bugfixes for calculating spectra and for misplaced/missing openMP statements

  • Property svn:keywords set to Id
File size: 66.9 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! Bugfix for misplaced and missing openMP directives in r2155
23!
24! Former revisions:
25! -----------------
26! $Id: prognostic_equations.f90 2192 2017-03-22 04:14:10Z 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!-- If required, calculate cloud microphysics
397    IF ( cloud_physics  .AND.  .NOT. microphysics_sat_adjust  .AND.            &
398         ( intermediate_timestep_count == 1  .OR.                              &
399           call_microphysics_at_all_substeps ) )                               &
400    THEN
401       !$OMP PARALLEL private (i,j)
402       !$OMP DO
403       DO  i = nxlg, nxrg
404          DO  j = nysg, nyng
405             CALL microphysics_control( i, j )
406           ENDDO
407       ENDDO
408       !$OMP END PARALLEL
409    ENDIF
410
411!
412!-- Loop over all prognostic equations
413    !$OMP PARALLEL private (i,i_omp_start,j,k,loop_start,tn)
414
415    !$  tn = omp_get_thread_num()
416    loop_start = .TRUE.
417
418    !$OMP DO
419    DO  i = nxl, nxr
420
421!
422!--    Store the first loop index. It differs for each thread and is required
423!--    later in advec_ws
424       IF ( loop_start )  THEN
425          loop_start  = .FALSE.
426          i_omp_start = i
427       ENDIF
428
429       DO  j = nys, nyn
430!
431!--       Tendency terms for u-velocity component
432          IF ( .NOT. outflow_l  .OR.  i > nxl )  THEN
433
434             tend(:,j,i) = 0.0_wp
435             IF ( timestep_scheme(1:5) == 'runge' )  THEN
436                IF ( ws_scheme_mom )  THEN
437                   CALL advec_u_ws( i, j, i_omp_start, tn )
438                ELSE
439                   CALL advec_u_pw( i, j )
440                ENDIF
441             ELSE
442                CALL advec_u_up( i, j )
443             ENDIF
444             CALL diffusion_u( i, j )
445             CALL coriolis( i, j, 1 )
446             IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
447                CALL buoyancy( i, j, pt, 1 )
448             ENDIF
449
450!
451!--          Drag by plant canopy
452             IF ( plant_canopy )  CALL pcm_tendency( i, j, 1 )
453
454!
455!--          External pressure gradient
456             IF ( dp_external )  THEN
457                DO  k = dp_level_ind_b+1, nzt
458                   tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
459                ENDDO
460             ENDIF
461
462!
463!--          Nudging
464             IF ( nudging )  CALL nudge( i, j, simulated_time, 'u' )
465
466!
467!--          Forces by wind turbines
468             IF ( wind_turbine )  CALL wtm_tendencies( i, j, 1 )
469
470             CALL user_actions( i, j, 'u-tendency' )
471!
472!--          Prognostic equation for u-velocity component
473             DO  k = nzb_u_inner(j,i)+1, nzt
474                u_p(k,j,i) = u(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +       &
475                                                  tsc(3) * tu_m(k,j,i) )       &
476                                      - tsc(5) * rdf(k) * ( u(k,j,i) - u_init(k) )
477             ENDDO
478
479!
480!--          Calculate tendencies for the next Runge-Kutta step
481             IF ( timestep_scheme(1:5) == 'runge' )  THEN
482                IF ( intermediate_timestep_count == 1 )  THEN
483                   DO  k = nzb_u_inner(j,i)+1, nzt
484                      tu_m(k,j,i) = tend(k,j,i)
485                   ENDDO
486                ELSEIF ( intermediate_timestep_count < &
487                         intermediate_timestep_count_max )  THEN
488                   DO  k = nzb_u_inner(j,i)+1, nzt
489                      tu_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tu_m(k,j,i)
490                   ENDDO
491                ENDIF
492             ENDIF
493
494          ENDIF
495
496!
497!--       Tendency terms for v-velocity component
498          IF ( .NOT. outflow_s  .OR.  j > nys )  THEN
499
500             tend(:,j,i) = 0.0_wp
501             IF ( timestep_scheme(1:5) == 'runge' )  THEN
502                IF ( ws_scheme_mom )  THEN
503                    CALL advec_v_ws( i, j, i_omp_start, tn )
504                ELSE
505                    CALL advec_v_pw( i, j )
506                ENDIF
507             ELSE
508                CALL advec_v_up( i, j )
509             ENDIF
510             CALL diffusion_v( i, j )
511             CALL coriolis( i, j, 2 )
512
513!
514!--          Drag by plant canopy
515             IF ( plant_canopy )  CALL pcm_tendency( i, j, 2 )
516
517!
518!--          External pressure gradient
519             IF ( dp_external )  THEN
520                DO  k = dp_level_ind_b+1, nzt
521                   tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
522                ENDDO
523             ENDIF
524
525!
526!--          Nudging
527             IF ( nudging )  CALL nudge( i, j, simulated_time, 'v' )
528
529!
530!--          Forces by wind turbines
531             IF ( wind_turbine )  CALL wtm_tendencies( i, j, 2 )
532
533             CALL user_actions( i, j, 'v-tendency' )
534!
535!--          Prognostic equation for v-velocity component
536             DO  k = nzb_v_inner(j,i)+1, nzt
537                v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +       &
538                                                  tsc(3) * tv_m(k,j,i) )       &
539                                      - tsc(5) * rdf(k) * ( v(k,j,i) - v_init(k) )
540             ENDDO
541
542!
543!--          Calculate tendencies for the next Runge-Kutta step
544             IF ( timestep_scheme(1:5) == 'runge' )  THEN
545                IF ( intermediate_timestep_count == 1 )  THEN
546                   DO  k = nzb_v_inner(j,i)+1, nzt
547                      tv_m(k,j,i) = tend(k,j,i)
548                   ENDDO
549                ELSEIF ( intermediate_timestep_count < &
550                         intermediate_timestep_count_max )  THEN
551                   DO  k = nzb_v_inner(j,i)+1, nzt
552                      tv_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tv_m(k,j,i)
553                   ENDDO
554                ENDIF
555             ENDIF
556
557          ENDIF
558
559!
560!--       Tendency terms for w-velocity component
561          tend(:,j,i) = 0.0_wp
562          IF ( timestep_scheme(1:5) == 'runge' )  THEN
563             IF ( ws_scheme_mom )  THEN
564                CALL advec_w_ws( i, j, i_omp_start, tn )
565             ELSE
566                CALL advec_w_pw( i, j )
567             END IF
568          ELSE
569             CALL advec_w_up( i, j )
570          ENDIF
571          CALL diffusion_w( i, j )
572          CALL coriolis( i, j, 3 )
573
574          IF ( .NOT. neutral )  THEN
575             IF ( ocean )  THEN
576                CALL buoyancy( i, j, rho_ocean, 3 )
577             ELSE
578                IF ( .NOT. humidity )  THEN
579                   CALL buoyancy( i, j, pt, 3 )
580                ELSE
581                   CALL buoyancy( i, j, vpt, 3 )
582                ENDIF
583             ENDIF
584          ENDIF
585
586!
587!--       Drag by plant canopy
588          IF ( plant_canopy )  CALL pcm_tendency( i, j, 3 )
589
590!
591!--       Forces by wind turbines
592          IF ( wind_turbine )  CALL wtm_tendencies( i, j, 3 )
593
594          CALL user_actions( i, j, 'w-tendency' )
595
596!
597!--       Prognostic equation for w-velocity component
598          DO  k = nzb_w_inner(j,i)+1, nzt-1
599             w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +          &
600                                               tsc(3) * tw_m(k,j,i) )          &
601                                   - tsc(5) * rdf(k) * w(k,j,i)
602          ENDDO
603
604!
605!--       Calculate tendencies for the next Runge-Kutta step
606          IF ( timestep_scheme(1:5) == 'runge' )  THEN
607             IF ( intermediate_timestep_count == 1 )  THEN
608                DO  k = nzb_w_inner(j,i)+1, nzt-1
609                   tw_m(k,j,i) = tend(k,j,i)
610                ENDDO
611             ELSEIF ( intermediate_timestep_count < &
612                      intermediate_timestep_count_max )  THEN
613                DO  k = nzb_w_inner(j,i)+1, nzt-1
614                   tw_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tw_m(k,j,i)
615                ENDDO
616             ENDIF
617          ENDIF
618
619!
620!--       If required, compute prognostic equation for potential temperature
621          IF ( .NOT. neutral )  THEN
622!
623!--          Tendency terms for potential temperature
624             tend(:,j,i) = 0.0_wp
625             IF ( timestep_scheme(1:5) == 'runge' )  THEN
626                   IF ( ws_scheme_sca )  THEN
627                      CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, diss_s_pt, &
628                                       flux_l_pt, diss_l_pt, i_omp_start, tn )
629                   ELSE
630                      CALL advec_s_pw( i, j, pt )
631                   ENDIF
632             ELSE
633                CALL advec_s_up( i, j, pt )
634             ENDIF
635             CALL diffusion_s( i, j, pt, shf, tswst, wall_heatflux )
636
637!
638!--          Tendency pt from wall heat flux from urban surface
639             IF ( urban_surface )  THEN
640                CALL usm_wall_heat_flux( i, j )
641             ENDIF
642
643!
644!--          If required compute heating/cooling due to long wave radiation
645!--          processes
646             IF ( cloud_top_radiation )  THEN
647                CALL calc_radiation( i, j )
648             ENDIF
649
650!
651!--          Consideration of heat sources within the plant canopy
652             IF ( plant_canopy  .AND.  cthf /= 0.0_wp )  THEN
653                CALL pcm_tendency( i, j, 4 )
654             ENDIF
655
656!
657!--          Large scale advection
658             IF ( large_scale_forcing )  THEN
659                CALL ls_advec( i, j, simulated_time, 'pt' )
660             ENDIF
661
662!
663!--          Nudging
664             IF ( nudging )  CALL nudge( i, j, simulated_time, 'pt' )
665
666!
667!--          If required, compute effect of large-scale subsidence/ascent
668             IF ( large_scale_subsidence  .AND.                                &
669                  .NOT. use_subsidence_tendencies )  THEN
670                CALL subsidence( i, j, tend, pt, pt_init, 2 )
671             ENDIF
672
673!
674!--          If required, add tendency due to radiative heating/cooling
675             IF ( radiation  .AND.                                             &
676                  simulated_time > skip_time_do_radiation )  THEN
677                CALL radiation_tendency ( i, j, tend )
678             ENDIF
679
680
681             CALL user_actions( i, j, 'pt-tendency' )
682
683!
684!--          Prognostic equation for potential temperature
685             DO  k = nzb_s_inner(j,i)+1, nzt
686                pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +     &
687                                                    tsc(3) * tpt_m(k,j,i) )    &
688                                        - tsc(5) * ( pt(k,j,i) - pt_init(k) ) *&
689                                          ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )
690             ENDDO
691
692!
693!--          Calculate tendencies for the next Runge-Kutta step
694             IF ( timestep_scheme(1:5) == 'runge' )  THEN
695                IF ( intermediate_timestep_count == 1 )  THEN
696                   DO  k = nzb_s_inner(j,i)+1, nzt
697                      tpt_m(k,j,i) = tend(k,j,i)
698                   ENDDO
699                ELSEIF ( intermediate_timestep_count < &
700                         intermediate_timestep_count_max )  THEN
701                   DO  k = nzb_s_inner(j,i)+1, nzt
702                      tpt_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
703                                      5.3125_wp * tpt_m(k,j,i)
704                   ENDDO
705                ENDIF
706             ENDIF
707
708          ENDIF
709
710!
711!--       If required, compute prognostic equation for salinity
712          IF ( ocean )  THEN
713
714!
715!--          Tendency-terms for salinity
716             tend(:,j,i) = 0.0_wp
717             IF ( timestep_scheme(1:5) == 'runge' ) &
718             THEN
719                IF ( ws_scheme_sca )  THEN
720                    CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa,  &
721                                diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn  )
722                ELSE
723                    CALL advec_s_pw( i, j, sa )
724                ENDIF
725             ELSE
726                CALL advec_s_up( i, j, sa )
727             ENDIF
728             CALL diffusion_s( i, j, sa, saswsb, saswst, wall_salinityflux )
729
730             CALL user_actions( i, j, 'sa-tendency' )
731
732!
733!--          Prognostic equation for salinity
734             DO  k = nzb_s_inner(j,i)+1, nzt
735                sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +     &
736                                                    tsc(3) * tsa_m(k,j,i) )    &
737                                        - tsc(5) * rdf_sc(k) *                 &
738                                          ( sa(k,j,i) - sa_init(k) )
739                IF ( sa_p(k,j,i) < 0.0_wp )  sa_p(k,j,i) = 0.1_wp * sa(k,j,i)
740             ENDDO
741
742!
743!--          Calculate tendencies for the next Runge-Kutta step
744             IF ( timestep_scheme(1:5) == 'runge' )  THEN
745                IF ( intermediate_timestep_count == 1 )  THEN
746                   DO  k = nzb_s_inner(j,i)+1, nzt
747                      tsa_m(k,j,i) = tend(k,j,i)
748                   ENDDO
749                ELSEIF ( intermediate_timestep_count < &
750                         intermediate_timestep_count_max )  THEN
751                   DO  k = nzb_s_inner(j,i)+1, nzt
752                      tsa_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
753                                      5.3125_wp * tsa_m(k,j,i)
754                   ENDDO
755                ENDIF
756             ENDIF
757
758!
759!--          Calculate density by the equation of state for seawater
760             CALL eqn_state_seawater( i, j )
761
762          ENDIF
763
764!
765!--       If required, compute prognostic equation for total water content
766          IF ( humidity )  THEN
767
768!
769!--          Tendency-terms for total water content / scalar
770             tend(:,j,i) = 0.0_wp
771             IF ( timestep_scheme(1:5) == 'runge' ) &
772             THEN
773                IF ( ws_scheme_sca )  THEN
774                   CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
775                                diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
776                ELSE
777                   CALL advec_s_pw( i, j, q )
778                ENDIF
779             ELSE
780                CALL advec_s_up( i, j, q )
781             ENDIF
782             CALL diffusion_s( i, j, q, qsws, qswst, wall_qflux )
783
784!
785!--          Sink or source of humidity due to canopy elements
786             IF ( plant_canopy )  CALL pcm_tendency( i, j, 5 )
787
788!
789!--          Large scale advection
790             IF ( large_scale_forcing )  THEN
791                CALL ls_advec( i, j, simulated_time, 'q' )
792             ENDIF
793
794!
795!--          Nudging
796             IF ( nudging )  CALL nudge( i, j, simulated_time, 'q' )
797
798!
799!--          If required compute influence of large-scale subsidence/ascent
800             IF ( large_scale_subsidence  .AND.                                &
801                  .NOT. use_subsidence_tendencies )  THEN
802                CALL subsidence( i, j, tend, q, q_init, 3 )
803             ENDIF
804
805             CALL user_actions( i, j, 'q-tendency' )
806
807!
808!--          Prognostic equation for total water content / scalar
809             DO  k = nzb_s_inner(j,i)+1, nzt
810                q_p(k,j,i) = q(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +       &
811                                                  tsc(3) * tq_m(k,j,i) )       &
812                                      - tsc(5) * rdf_sc(k) *                   &
813                                        ( q(k,j,i) - q_init(k) )
814                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
815             ENDDO
816
817!
818!--          Calculate tendencies for the next Runge-Kutta step
819             IF ( timestep_scheme(1:5) == 'runge' )  THEN
820                IF ( intermediate_timestep_count == 1 )  THEN
821                   DO  k = nzb_s_inner(j,i)+1, nzt
822                      tq_m(k,j,i) = tend(k,j,i)
823                   ENDDO
824                ELSEIF ( intermediate_timestep_count < &
825                         intermediate_timestep_count_max )  THEN
826                   DO  k = nzb_s_inner(j,i)+1, nzt
827                      tq_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
828                                     5.3125_wp * tq_m(k,j,i)
829                   ENDDO
830                ENDIF
831             ENDIF
832
833!
834!--          If required, calculate prognostic equations for rain water content
835!--          and rain drop concentration
836             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
837!
838!--             Calculate prognostic equation for rain water content
839                tend(:,j,i) = 0.0_wp
840                IF ( timestep_scheme(1:5) == 'runge' ) &
841                THEN
842                   IF ( ws_scheme_sca )  THEN
843                      CALL advec_s_ws( i, j, qr, 'qr', flux_s_qr,       &
844                                       diss_s_qr, flux_l_qr, diss_l_qr, &
845                                       i_omp_start, tn )
846                   ELSE
847                      CALL advec_s_pw( i, j, qr )
848                   ENDIF
849                ELSE
850                   CALL advec_s_up( i, j, qr )
851                ENDIF
852                CALL diffusion_s( i, j, qr, qrsws, qrswst, wall_qrflux )
853
854!
855!--             Prognostic equation for rain water content
856                DO  k = nzb_s_inner(j,i)+1, nzt
857                   qr_p(k,j,i) = qr(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +  &
858                                                       tsc(3) * tqr_m(k,j,i) ) &
859                                           - tsc(5) * rdf_sc(k) * qr(k,j,i)
860                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
861                ENDDO
862!
863!--             Calculate tendencies for the next Runge-Kutta step
864                IF ( timestep_scheme(1:5) == 'runge' )  THEN
865                   IF ( intermediate_timestep_count == 1 )  THEN
866                      DO  k = nzb_s_inner(j,i)+1, nzt
867                         tqr_m(k,j,i) = tend(k,j,i)
868                      ENDDO
869                   ELSEIF ( intermediate_timestep_count < &
870                            intermediate_timestep_count_max )  THEN
871                      DO  k = nzb_s_inner(j,i)+1, nzt
872                         tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
873                                        5.3125_wp * tqr_m(k,j,i)
874                      ENDDO
875                   ENDIF
876                ENDIF
877
878!
879!--             Calculate prognostic equation for rain drop concentration.
880                tend(:,j,i) = 0.0_wp
881                IF ( timestep_scheme(1:5) == 'runge' )  THEN
882                   IF ( ws_scheme_sca )  THEN
883                      CALL advec_s_ws( i, j, nr, 'nr', flux_s_nr,    &
884                                    diss_s_nr, flux_l_nr, diss_l_nr, &
885                                    i_omp_start, tn )
886                   ELSE
887                      CALL advec_s_pw( i, j, nr )
888                   ENDIF
889                ELSE
890                   CALL advec_s_up( i, j, nr )
891                ENDIF
892                CALL diffusion_s( i, j, nr, nrsws, nrswst, wall_nrflux )
893
894!
895!--             Prognostic equation for rain drop concentration
896                DO  k = nzb_s_inner(j,i)+1, nzt
897                   nr_p(k,j,i) = nr(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +  &
898                                                       tsc(3) * tnr_m(k,j,i) ) &
899                                           - tsc(5) * rdf_sc(k) * nr(k,j,i)
900                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
901                ENDDO
902!
903!--             Calculate tendencies for the next Runge-Kutta step
904                IF ( timestep_scheme(1:5) == 'runge' )  THEN
905                   IF ( intermediate_timestep_count == 1 )  THEN
906                      DO  k = nzb_s_inner(j,i)+1, nzt
907                         tnr_m(k,j,i) = tend(k,j,i)
908                      ENDDO
909                   ELSEIF ( intermediate_timestep_count < &
910                            intermediate_timestep_count_max )  THEN
911                      DO  k = nzb_s_inner(j,i)+1, nzt
912                         tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
913                                        5.3125_wp * tnr_m(k,j,i)
914                      ENDDO
915                   ENDIF
916                ENDIF
917
918             ENDIF
919
920          ENDIF
921
922!
923!--       If required, compute prognostic equation for scalar
924          IF ( passive_scalar )  THEN
925!
926!--          Tendency-terms for total water content / scalar
927             tend(:,j,i) = 0.0_wp
928             IF ( timestep_scheme(1:5) == 'runge' ) &
929             THEN
930                IF ( ws_scheme_sca )  THEN
931                   CALL advec_s_ws( i, j, s, 's', flux_s_s, &
932                                diss_s_s, flux_l_s, diss_l_s, i_omp_start, tn )
933                ELSE
934                   CALL advec_s_pw( i, j, s )
935                ENDIF
936             ELSE
937                CALL advec_s_up( i, j, s )
938             ENDIF
939             CALL diffusion_s( i, j, s, ssws, sswst, wall_sflux )
940
941!
942!--          Sink or source of scalar concentration due to canopy elements
943             IF ( plant_canopy )  CALL pcm_tendency( i, j, 7 )
944
945!
946!--          Large scale advection, still need to be extended for scalars
947!              IF ( large_scale_forcing )  THEN
948!                 CALL ls_advec( i, j, simulated_time, 's' )
949!              ENDIF
950
951!
952!--          Nudging, still need to be extended for scalars
953!              IF ( nudging )  CALL nudge( i, j, simulated_time, 's' )
954
955!
956!--          If required compute influence of large-scale subsidence/ascent.
957!--          Note, the last argument is of no meaning in this case, as it is
958!--          only used in conjunction with large_scale_forcing, which is to
959!--          date not implemented for scalars.
960             IF ( large_scale_subsidence  .AND.                                &
961                  .NOT. use_subsidence_tendencies  .AND.                       &
962                  .NOT. large_scale_forcing )  THEN
963                CALL subsidence( i, j, tend, s, s_init, 3 )
964             ENDIF
965
966             CALL user_actions( i, j, 's-tendency' )
967
968!
969!--          Prognostic equation for scalar
970             DO  k = nzb_s_inner(j,i)+1, nzt
971                s_p(k,j,i) = s(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +       &
972                                                  tsc(3) * ts_m(k,j,i) )       &
973                                      - tsc(5) * rdf_sc(k) *                   &
974                                        ( s(k,j,i) - s_init(k) )
975                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
976             ENDDO
977
978!
979!--          Calculate tendencies for the next Runge-Kutta step
980             IF ( timestep_scheme(1:5) == 'runge' )  THEN
981                IF ( intermediate_timestep_count == 1 )  THEN
982                   DO  k = nzb_s_inner(j,i)+1, nzt
983                      ts_m(k,j,i) = tend(k,j,i)
984                   ENDDO
985                ELSEIF ( intermediate_timestep_count < &
986                         intermediate_timestep_count_max )  THEN
987                   DO  k = nzb_s_inner(j,i)+1, nzt
988                      ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
989                                     5.3125_wp * ts_m(k,j,i)
990                   ENDDO
991                ENDIF
992             ENDIF
993
994          ENDIF
995!
996!--       If required, compute prognostic equation for turbulent kinetic
997!--       energy (TKE)
998          IF ( .NOT. constant_diffusion )  THEN
999
1000!
1001!--          Tendency-terms for TKE
1002             tend(:,j,i) = 0.0_wp
1003             IF ( timestep_scheme(1:5) == 'runge'  &
1004                 .AND.  .NOT. use_upstream_for_tke )  THEN
1005                 IF ( ws_scheme_sca )  THEN
1006                     CALL advec_s_ws( i, j, e, 'e', flux_s_e, diss_s_e, &
1007                                      flux_l_e, diss_l_e , i_omp_start, tn )
1008                 ELSE
1009                     CALL advec_s_pw( i, j, e )
1010                 ENDIF
1011             ELSE
1012                CALL advec_s_up( i, j, e )
1013             ENDIF
1014             IF ( .NOT. humidity )  THEN
1015                IF ( ocean )  THEN
1016                   CALL diffusion_e( i, j, prho, prho_reference )
1017                ELSE
1018                   CALL diffusion_e( i, j, pt, pt_reference )
1019                ENDIF
1020             ELSE
1021                CALL diffusion_e( i, j, vpt, pt_reference )
1022             ENDIF
1023             CALL production_e( i, j )
1024
1025!
1026!--          Additional sink term for flows through plant canopies
1027             IF ( plant_canopy )  CALL pcm_tendency( i, j, 6 )
1028
1029             CALL user_actions( i, j, 'e-tendency' )
1030
1031!
1032!--          Prognostic equation for TKE.
1033!--          Eliminate negative TKE values, which can occur due to numerical
1034!--          reasons in the course of the integration. In such cases the old
1035!--          TKE value is reduced by 90%.
1036             DO  k = nzb_s_inner(j,i)+1, nzt
1037                e_p(k,j,i) = e(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +       &
1038                                                  tsc(3) * te_m(k,j,i) )
1039                IF ( e_p(k,j,i) < 0.0_wp )  e_p(k,j,i) = 0.1_wp * e(k,j,i)
1040             ENDDO
1041
1042!
1043!--          Calculate tendencies for the next Runge-Kutta step
1044             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1045                IF ( intermediate_timestep_count == 1 )  THEN
1046                   DO  k = nzb_s_inner(j,i)+1, nzt
1047                      te_m(k,j,i) = tend(k,j,i)
1048                   ENDDO
1049                ELSEIF ( intermediate_timestep_count < &
1050                         intermediate_timestep_count_max )  THEN
1051                   DO  k = nzb_s_inner(j,i)+1, nzt
1052                      te_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1053                                     5.3125_wp * te_m(k,j,i)
1054                   ENDDO
1055                ENDIF
1056             ENDIF
1057
1058          ENDIF   ! TKE equation
1059
1060       ENDDO
1061    ENDDO
1062    !$OMP END PARALLEL
1063
1064    CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1065
1066
1067 END SUBROUTINE prognostic_equations_cache
1068
1069
1070!------------------------------------------------------------------------------!
1071! Description:
1072! ------------
1073!> Version for vector machines
1074!------------------------------------------------------------------------------!
1075
1076 SUBROUTINE prognostic_equations_vector
1077
1078
1079    IMPLICIT NONE
1080
1081    INTEGER(iwp) ::  i    !<
1082    INTEGER(iwp) ::  j    !<
1083    INTEGER(iwp) ::  k    !<
1084
1085    REAL(wp)     ::  sbt  !<
1086
1087
1088!
1089!-- If required, calculate cloud microphysical impacts
1090    IF ( cloud_physics  .AND.  .NOT. microphysics_sat_adjust  .AND.            &
1091         ( intermediate_timestep_count == 1  .OR.                              &
1092           call_microphysics_at_all_substeps )                                 &
1093       )  THEN
1094       CALL cpu_log( log_point(51), 'microphysics', 'start' )
1095       CALL microphysics_control
1096       CALL cpu_log( log_point(51), 'microphysics', 'stop' )
1097    ENDIF
1098
1099!
1100!-- u-velocity component
1101    CALL cpu_log( log_point(5), 'u-equation', 'start' )
1102
1103    tend = 0.0_wp
1104    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1105       IF ( ws_scheme_mom )  THEN
1106          CALL advec_u_ws
1107       ELSE
1108          CALL advec_u_pw
1109       ENDIF
1110    ELSE
1111       CALL advec_u_up
1112    ENDIF
1113    CALL diffusion_u
1114    CALL coriolis( 1 )
1115    IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
1116       CALL buoyancy( pt, 1 )
1117    ENDIF
1118
1119!
1120!-- Drag by plant canopy
1121    IF ( plant_canopy )  CALL pcm_tendency( 1 )
1122
1123!
1124!-- External pressure gradient
1125    IF ( dp_external )  THEN
1126       DO  i = nxlu, nxr
1127          DO  j = nys, nyn
1128             DO  k = dp_level_ind_b+1, nzt
1129                tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
1130             ENDDO
1131          ENDDO
1132       ENDDO
1133    ENDIF
1134
1135!
1136!-- Nudging
1137    IF ( nudging )  CALL nudge( simulated_time, 'u' )
1138
1139!
1140!-- Forces by wind turbines
1141    IF ( wind_turbine )  CALL wtm_tendencies( 1 )
1142
1143    CALL user_actions( 'u-tendency' )
1144
1145!
1146!-- Prognostic equation for u-velocity component
1147    DO  i = nxlu, nxr
1148       DO  j = nys, nyn
1149          DO  k = nzb_u_inner(j,i)+1, nzt
1150             u_p(k,j,i) = u(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +          &
1151                                               tsc(3) * tu_m(k,j,i) )          &
1152                                   - tsc(5) * rdf(k) * ( u(k,j,i) - u_init(k) )
1153          ENDDO
1154       ENDDO
1155    ENDDO
1156
1157!
1158!-- Calculate tendencies for the next Runge-Kutta step
1159    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1160       IF ( intermediate_timestep_count == 1 )  THEN
1161          DO  i = nxlu, nxr
1162             DO  j = nys, nyn
1163                DO  k = nzb_u_inner(j,i)+1, nzt
1164                   tu_m(k,j,i) = tend(k,j,i)
1165                ENDDO
1166             ENDDO
1167          ENDDO
1168       ELSEIF ( intermediate_timestep_count < &
1169                intermediate_timestep_count_max )  THEN
1170          DO  i = nxlu, nxr
1171             DO  j = nys, nyn
1172                DO  k = nzb_u_inner(j,i)+1, nzt
1173                   tu_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tu_m(k,j,i)
1174                ENDDO
1175             ENDDO
1176          ENDDO
1177       ENDIF
1178    ENDIF
1179
1180    CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1181
1182!
1183!-- v-velocity component
1184    CALL cpu_log( log_point(6), 'v-equation', 'start' )
1185
1186    tend = 0.0_wp
1187    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1188       IF ( ws_scheme_mom )  THEN
1189          CALL advec_v_ws
1190       ELSE
1191          CALL advec_v_pw
1192       END IF
1193    ELSE
1194       CALL advec_v_up
1195    ENDIF
1196    CALL diffusion_v
1197    CALL coriolis( 2 )
1198
1199!
1200!-- Drag by plant canopy
1201    IF ( plant_canopy )  CALL pcm_tendency( 2 )
1202
1203!
1204!-- External pressure gradient
1205    IF ( dp_external )  THEN
1206       DO  i = nxl, nxr
1207          DO  j = nysv, nyn
1208             DO  k = dp_level_ind_b+1, nzt
1209                tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1210             ENDDO
1211          ENDDO
1212       ENDDO
1213    ENDIF
1214
1215!
1216!-- Nudging
1217    IF ( nudging )  CALL nudge( simulated_time, 'v' )
1218
1219!
1220!-- Forces by wind turbines
1221    IF ( wind_turbine )  CALL wtm_tendencies( 2 )
1222
1223    CALL user_actions( 'v-tendency' )
1224
1225!
1226!-- Prognostic equation for v-velocity component
1227    DO  i = nxl, nxr
1228       DO  j = nysv, nyn
1229          DO  k = nzb_v_inner(j,i)+1, nzt
1230             v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +          &
1231                                               tsc(3) * tv_m(k,j,i) )          &
1232                                   - tsc(5) * rdf(k) * ( v(k,j,i) - v_init(k) )
1233          ENDDO
1234       ENDDO
1235    ENDDO
1236
1237!
1238!-- Calculate tendencies for the next Runge-Kutta step
1239    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1240       IF ( intermediate_timestep_count == 1 )  THEN
1241          DO  i = nxl, nxr
1242             DO  j = nysv, nyn
1243                DO  k = nzb_v_inner(j,i)+1, nzt
1244                   tv_m(k,j,i) = tend(k,j,i)
1245                ENDDO
1246             ENDDO
1247          ENDDO
1248       ELSEIF ( intermediate_timestep_count < &
1249                intermediate_timestep_count_max )  THEN
1250          DO  i = nxl, nxr
1251             DO  j = nysv, nyn
1252                DO  k = nzb_v_inner(j,i)+1, nzt
1253                   tv_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tv_m(k,j,i)
1254                ENDDO
1255             ENDDO
1256          ENDDO
1257       ENDIF
1258    ENDIF
1259
1260    CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1261
1262!
1263!-- w-velocity component
1264    CALL cpu_log( log_point(7), 'w-equation', 'start' )
1265
1266    tend = 0.0_wp
1267    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1268       IF ( ws_scheme_mom )  THEN
1269          CALL advec_w_ws
1270       ELSE
1271          CALL advec_w_pw
1272       ENDIF
1273    ELSE
1274       CALL advec_w_up
1275    ENDIF
1276    CALL diffusion_w
1277    CALL coriolis( 3 )
1278
1279    IF ( .NOT. neutral )  THEN
1280       IF ( ocean )  THEN
1281          CALL buoyancy( rho_ocean, 3 )
1282       ELSE
1283          IF ( .NOT. humidity )  THEN
1284             CALL buoyancy( pt, 3 )
1285          ELSE
1286             CALL buoyancy( vpt, 3 )
1287          ENDIF
1288       ENDIF
1289    ENDIF
1290
1291!
1292!-- Drag by plant canopy
1293    IF ( plant_canopy )  CALL pcm_tendency( 3 )
1294
1295!
1296!-- Forces by wind turbines
1297    IF ( wind_turbine )  CALL wtm_tendencies( 3 )
1298
1299    CALL user_actions( 'w-tendency' )
1300
1301!
1302!-- Prognostic equation for w-velocity component
1303    DO  i = nxl, nxr
1304       DO  j = nys, nyn
1305          DO  k = nzb_w_inner(j,i)+1, nzt-1
1306             w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +          &
1307                                               tsc(3) * tw_m(k,j,i) )          &
1308                                   - tsc(5) * rdf(k) * w(k,j,i)
1309          ENDDO
1310       ENDDO
1311    ENDDO
1312
1313!
1314!-- Calculate tendencies for the next Runge-Kutta step
1315    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1316       IF ( intermediate_timestep_count == 1 )  THEN
1317          DO  i = nxl, nxr
1318             DO  j = nys, nyn
1319                DO  k = nzb_w_inner(j,i)+1, nzt-1
1320                   tw_m(k,j,i) = tend(k,j,i)
1321                ENDDO
1322             ENDDO
1323          ENDDO
1324       ELSEIF ( intermediate_timestep_count < &
1325                intermediate_timestep_count_max )  THEN
1326          DO  i = nxl, nxr
1327             DO  j = nys, nyn
1328                DO  k = nzb_w_inner(j,i)+1, nzt-1
1329                   tw_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tw_m(k,j,i)
1330                ENDDO
1331             ENDDO
1332          ENDDO
1333       ENDIF
1334    ENDIF
1335
1336    CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1337
1338
1339!
1340!-- If required, compute prognostic equation for potential temperature
1341    IF ( .NOT. neutral )  THEN
1342
1343       CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1344
1345!
1346!--    pt-tendency terms with communication
1347       sbt = tsc(2)
1348       IF ( scalar_advec == 'bc-scheme' )  THEN
1349
1350          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1351!
1352!--          Bott-Chlond scheme always uses Euler time step. Thus:
1353             sbt = 1.0_wp
1354          ENDIF
1355          tend = 0.0_wp
1356          CALL advec_s_bc( pt, 'pt' )
1357
1358       ENDIF
1359
1360!
1361!--    pt-tendency terms with no communication
1362       IF ( scalar_advec /= 'bc-scheme' )  THEN
1363          tend = 0.0_wp
1364          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1365             IF ( ws_scheme_sca )  THEN
1366                CALL advec_s_ws( pt, 'pt' )
1367             ELSE
1368                CALL advec_s_pw( pt )
1369             ENDIF
1370          ELSE
1371             CALL advec_s_up( pt )
1372          ENDIF
1373       ENDIF
1374
1375       CALL diffusion_s( pt, shf, tswst, wall_heatflux )
1376
1377!
1378!--    Tendency pt from wall heat flux from urban surface
1379       IF ( urban_surface )  THEN
1380          CALL usm_wall_heat_flux
1381       ENDIF
1382
1383!
1384!--    If required compute heating/cooling due to long wave radiation processes
1385       IF ( cloud_top_radiation )  THEN
1386          CALL calc_radiation
1387       ENDIF
1388
1389!
1390!--    Consideration of heat sources within the plant canopy
1391       IF ( plant_canopy .AND. ( cthf /= 0.0_wp ) ) THEN
1392          CALL pcm_tendency( 4 )
1393       ENDIF
1394
1395!
1396!--    Large scale advection
1397       IF ( large_scale_forcing )  THEN
1398          CALL ls_advec( simulated_time, 'pt' )
1399       ENDIF
1400
1401!
1402!--    Nudging
1403       IF ( nudging )  CALL nudge( simulated_time, 'pt' )
1404
1405!
1406!--    If required compute influence of large-scale subsidence/ascent
1407       IF ( large_scale_subsidence  .AND.                                      &
1408            .NOT. use_subsidence_tendencies )  THEN
1409          CALL subsidence( tend, pt, pt_init, 2 )
1410       ENDIF
1411
1412!
1413!--    If required, add tendency due to radiative heating/cooling
1414       IF ( radiation  .AND.                                                   &
1415            simulated_time > skip_time_do_radiation )  THEN
1416            CALL radiation_tendency ( tend )
1417       ENDIF
1418
1419       CALL user_actions( 'pt-tendency' )
1420
1421!
1422!--    Prognostic equation for potential temperature
1423       DO  i = nxl, nxr
1424          DO  j = nys, nyn
1425             DO  k = nzb_s_inner(j,i)+1, nzt
1426                pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +        &
1427                                                    tsc(3) * tpt_m(k,j,i) )    &
1428                                        - tsc(5) * ( pt(k,j,i) - pt_init(k) ) *&
1429                                          ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )
1430             ENDDO
1431          ENDDO
1432       ENDDO
1433
1434!
1435!--    Calculate tendencies for the next Runge-Kutta step
1436       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1437          IF ( intermediate_timestep_count == 1 )  THEN
1438             DO  i = nxl, nxr
1439                DO  j = nys, nyn
1440                   DO  k = nzb_s_inner(j,i)+1, nzt
1441                      tpt_m(k,j,i) = tend(k,j,i)
1442                   ENDDO
1443                ENDDO
1444             ENDDO
1445          ELSEIF ( intermediate_timestep_count < &
1446                   intermediate_timestep_count_max )  THEN
1447             DO  i = nxl, nxr
1448                DO  j = nys, nyn
1449                   DO  k = nzb_s_inner(j,i)+1, nzt
1450                      tpt_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1451                                      5.3125_wp * tpt_m(k,j,i)
1452                   ENDDO
1453                ENDDO
1454             ENDDO
1455          ENDIF
1456       ENDIF
1457
1458       CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1459
1460    ENDIF
1461
1462!
1463!-- If required, compute prognostic equation for salinity
1464    IF ( ocean )  THEN
1465
1466       CALL cpu_log( log_point(37), 'sa-equation', 'start' )
1467
1468!
1469!--    sa-tendency terms with communication
1470       sbt = tsc(2)
1471       IF ( scalar_advec == 'bc-scheme' )  THEN
1472
1473          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1474!
1475!--          Bott-Chlond scheme always uses Euler time step. Thus:
1476             sbt = 1.0_wp
1477          ENDIF
1478          tend = 0.0_wp
1479          CALL advec_s_bc( sa, 'sa' )
1480
1481       ENDIF
1482
1483!
1484!--    sa-tendency terms with no communication
1485       IF ( scalar_advec /= 'bc-scheme' )  THEN
1486          tend = 0.0_wp
1487          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1488             IF ( ws_scheme_sca )  THEN
1489                 CALL advec_s_ws( sa, 'sa' )
1490             ELSE
1491                 CALL advec_s_pw( sa )
1492             ENDIF
1493          ELSE
1494             CALL advec_s_up( sa )
1495          ENDIF
1496       ENDIF
1497
1498       CALL diffusion_s( sa, saswsb, saswst, wall_salinityflux )
1499
1500       CALL user_actions( 'sa-tendency' )
1501
1502!
1503!--    Prognostic equation for salinity
1504       DO  i = nxl, nxr
1505          DO  j = nys, nyn
1506             DO  k = nzb_s_inner(j,i)+1, nzt
1507                sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +        &
1508                                                    tsc(3) * tsa_m(k,j,i) )    &
1509                                        - tsc(5) * rdf_sc(k) *                 &
1510                                          ( sa(k,j,i) - sa_init(k) )
1511                IF ( sa_p(k,j,i) < 0.0_wp )  sa_p(k,j,i) = 0.1_wp * sa(k,j,i)
1512             ENDDO
1513          ENDDO
1514       ENDDO
1515
1516!
1517!--    Calculate tendencies for the next Runge-Kutta step
1518       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1519          IF ( intermediate_timestep_count == 1 )  THEN
1520             DO  i = nxl, nxr
1521                DO  j = nys, nyn
1522                   DO  k = nzb_s_inner(j,i)+1, nzt
1523                      tsa_m(k,j,i) = tend(k,j,i)
1524                   ENDDO
1525                ENDDO
1526             ENDDO
1527          ELSEIF ( intermediate_timestep_count < &
1528                   intermediate_timestep_count_max )  THEN
1529             DO  i = nxl, nxr
1530                DO  j = nys, nyn
1531                   DO  k = nzb_s_inner(j,i)+1, nzt
1532                      tsa_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1533                                      5.3125_wp * tsa_m(k,j,i)
1534                   ENDDO
1535                ENDDO
1536             ENDDO
1537          ENDIF
1538       ENDIF
1539
1540       CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
1541
1542!
1543!--    Calculate density by the equation of state for seawater
1544       CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
1545       CALL eqn_state_seawater
1546       CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
1547
1548    ENDIF
1549
1550!
1551!-- If required, compute prognostic equation for total water content
1552    IF ( humidity )  THEN
1553
1554       CALL cpu_log( log_point(29), 'q-equation', 'start' )
1555
1556!
1557!--    Scalar/q-tendency terms with communication
1558       sbt = tsc(2)
1559       IF ( scalar_advec == 'bc-scheme' )  THEN
1560
1561          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1562!
1563!--          Bott-Chlond scheme always uses Euler time step. Thus:
1564             sbt = 1.0_wp
1565          ENDIF
1566          tend = 0.0_wp
1567          CALL advec_s_bc( q, 'q' )
1568
1569       ENDIF
1570
1571!
1572!--    Scalar/q-tendency terms with no communication
1573       IF ( scalar_advec /= 'bc-scheme' )  THEN
1574          tend = 0.0_wp
1575          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1576             IF ( ws_scheme_sca )  THEN
1577                CALL advec_s_ws( q, 'q' )
1578             ELSE
1579                CALL advec_s_pw( q )
1580             ENDIF
1581          ELSE
1582             CALL advec_s_up( q )
1583          ENDIF
1584       ENDIF
1585
1586       CALL diffusion_s( q, qsws, qswst, wall_qflux )
1587
1588!
1589!--    Sink or source of humidity due to canopy elements
1590       IF ( plant_canopy ) CALL pcm_tendency( 5 )
1591
1592!
1593!--    Large scale advection
1594       IF ( large_scale_forcing )  THEN
1595          CALL ls_advec( simulated_time, 'q' )
1596       ENDIF
1597
1598!
1599!--    Nudging
1600       IF ( nudging )  CALL nudge( simulated_time, 'q' )
1601
1602!
1603!--    If required compute influence of large-scale subsidence/ascent
1604       IF ( large_scale_subsidence  .AND.                                      &
1605            .NOT. use_subsidence_tendencies )  THEN
1606         CALL subsidence( tend, q, q_init, 3 )
1607       ENDIF
1608
1609       CALL user_actions( 'q-tendency' )
1610
1611!
1612!--    Prognostic equation for total water content
1613       DO  i = nxl, nxr
1614          DO  j = nys, nyn
1615             DO  k = nzb_s_inner(j,i)+1, nzt
1616                q_p(k,j,i) = q(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +          &
1617                                                  tsc(3) * tq_m(k,j,i) )       &
1618                                      - tsc(5) * rdf_sc(k) *                   &
1619                                        ( q(k,j,i) - q_init(k) )
1620                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
1621             ENDDO
1622          ENDDO
1623       ENDDO
1624
1625!
1626!--    Calculate tendencies for the next Runge-Kutta step
1627       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1628          IF ( intermediate_timestep_count == 1 )  THEN
1629             DO  i = nxl, nxr
1630                DO  j = nys, nyn
1631                   DO  k = nzb_s_inner(j,i)+1, nzt
1632                      tq_m(k,j,i) = tend(k,j,i)
1633                   ENDDO
1634                ENDDO
1635             ENDDO
1636          ELSEIF ( intermediate_timestep_count < &
1637                   intermediate_timestep_count_max )  THEN
1638             DO  i = nxl, nxr
1639                DO  j = nys, nyn
1640                   DO  k = nzb_s_inner(j,i)+1, nzt
1641                      tq_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tq_m(k,j,i)
1642                   ENDDO
1643                ENDDO
1644             ENDDO
1645          ENDIF
1646       ENDIF
1647
1648       CALL cpu_log( log_point(29), 'q-equation', 'stop' )
1649
1650!
1651!--    If required, calculate prognostic equations for rain water content
1652!--    and rain drop concentration
1653       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
1654
1655          CALL cpu_log( log_point(52), 'qr-equation', 'start' )
1656
1657!
1658!--       Calculate prognostic equation for rain water content
1659          sbt = tsc(2)
1660          IF ( scalar_advec == 'bc-scheme' )  THEN
1661
1662             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1663!
1664!--             Bott-Chlond scheme always uses Euler time step. Thus:
1665                sbt = 1.0_wp
1666             ENDIF
1667             tend = 0.0_wp
1668             CALL advec_s_bc( qr, 'qr' )
1669
1670          ENDIF
1671
1672!
1673!--       qr-tendency terms with no communication
1674          IF ( scalar_advec /= 'bc-scheme' )  THEN
1675             tend = 0.0_wp
1676             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1677                IF ( ws_scheme_sca )  THEN
1678                   CALL advec_s_ws( qr, 'qr' )
1679                ELSE
1680                   CALL advec_s_pw( qr )
1681                ENDIF
1682             ELSE
1683                CALL advec_s_up( qr )
1684             ENDIF
1685          ENDIF
1686
1687          CALL diffusion_s( qr, qrsws, qrswst, wall_qrflux )
1688
1689!
1690!--       Prognostic equation for rain water content
1691          DO  i = nxl, nxr
1692             DO  j = nys, nyn
1693                DO  k = nzb_s_inner(j,i)+1, nzt
1694                   qr_p(k,j,i) = qr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +     &
1695                                                     tsc(3) * tqr_m(k,j,i) )   &
1696                                           - tsc(5) * rdf_sc(k) * qr(k,j,i)
1697                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
1698                ENDDO
1699             ENDDO
1700          ENDDO
1701
1702!
1703!--       Calculate tendencies for the next Runge-Kutta step
1704          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1705             IF ( intermediate_timestep_count == 1 )  THEN
1706                DO  i = nxl, nxr
1707                   DO  j = nys, nyn
1708                      DO  k = nzb_s_inner(j,i)+1, nzt
1709                         tqr_m(k,j,i) = tend(k,j,i)
1710                      ENDDO
1711                   ENDDO
1712                ENDDO
1713             ELSEIF ( intermediate_timestep_count < &
1714                      intermediate_timestep_count_max )  THEN
1715                DO  i = nxl, nxr
1716                   DO  j = nys, nyn
1717                      DO  k = nzb_s_inner(j,i)+1, nzt
1718                         tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * &
1719                                                                   tqr_m(k,j,i)
1720                      ENDDO
1721                   ENDDO
1722                ENDDO
1723             ENDIF
1724          ENDIF
1725
1726          CALL cpu_log( log_point(52), 'qr-equation', 'stop' )
1727          CALL cpu_log( log_point(53), 'nr-equation', 'start' )
1728
1729!
1730!--       Calculate prognostic equation for rain drop concentration
1731          sbt = tsc(2)
1732          IF ( scalar_advec == 'bc-scheme' )  THEN
1733
1734             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1735!
1736!--             Bott-Chlond scheme always uses Euler time step. Thus:
1737                sbt = 1.0_wp
1738             ENDIF
1739             tend = 0.0_wp
1740             CALL advec_s_bc( nr, 'nr' )
1741
1742          ENDIF
1743
1744!
1745!--       nr-tendency terms with no communication
1746          IF ( scalar_advec /= 'bc-scheme' )  THEN
1747             tend = 0.0_wp
1748             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1749                IF ( ws_scheme_sca )  THEN
1750                   CALL advec_s_ws( nr, 'nr' )
1751                ELSE
1752                   CALL advec_s_pw( nr )
1753                ENDIF
1754             ELSE
1755                CALL advec_s_up( nr )
1756             ENDIF
1757          ENDIF
1758
1759          CALL diffusion_s( nr, nrsws, nrswst, wall_nrflux )
1760
1761!
1762!--       Prognostic equation for rain drop concentration
1763          DO  i = nxl, nxr
1764             DO  j = nys, nyn
1765                DO  k = nzb_s_inner(j,i)+1, nzt
1766                   nr_p(k,j,i) = nr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +     &
1767                                                     tsc(3) * tnr_m(k,j,i) )   &
1768                                           - tsc(5) * rdf_sc(k) * nr(k,j,i)
1769                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
1770                ENDDO
1771             ENDDO
1772          ENDDO
1773
1774!
1775!--       Calculate tendencies for the next Runge-Kutta step
1776          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1777             IF ( intermediate_timestep_count == 1 )  THEN
1778                DO  i = nxl, nxr
1779                   DO  j = nys, nyn
1780                      DO  k = nzb_s_inner(j,i)+1, nzt
1781                         tnr_m(k,j,i) = tend(k,j,i)
1782                      ENDDO
1783                   ENDDO
1784                ENDDO
1785             ELSEIF ( intermediate_timestep_count < &
1786                      intermediate_timestep_count_max )  THEN
1787                DO  i = nxl, nxr
1788                   DO  j = nys, nyn
1789                      DO  k = nzb_s_inner(j,i)+1, nzt
1790                         tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * &
1791                                                                   tnr_m(k,j,i)
1792                      ENDDO
1793                   ENDDO
1794                ENDDO
1795             ENDIF
1796          ENDIF
1797
1798          CALL cpu_log( log_point(53), 'nr-equation', 'stop' )
1799
1800       ENDIF
1801
1802    ENDIF
1803!
1804!-- If required, compute prognostic equation for scalar
1805    IF ( passive_scalar )  THEN
1806
1807       CALL cpu_log( log_point(66), 's-equation', 'start' )
1808
1809!
1810!--    Scalar/q-tendency terms with communication
1811       sbt = tsc(2)
1812       IF ( scalar_advec == 'bc-scheme' )  THEN
1813
1814          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1815!
1816!--          Bott-Chlond scheme always uses Euler time step. Thus:
1817             sbt = 1.0_wp
1818          ENDIF
1819          tend = 0.0_wp
1820          CALL advec_s_bc( s, 's' )
1821
1822       ENDIF
1823
1824!
1825!--    Scalar/q-tendency terms with no communication
1826       IF ( scalar_advec /= 'bc-scheme' )  THEN
1827          tend = 0.0_wp
1828          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1829             IF ( ws_scheme_sca )  THEN
1830                CALL advec_s_ws( s, 's' )
1831             ELSE
1832                CALL advec_s_pw( s )
1833             ENDIF
1834          ELSE
1835             CALL advec_s_up( s )
1836          ENDIF
1837       ENDIF
1838
1839       CALL diffusion_s( s, ssws, sswst, wall_sflux )
1840
1841!
1842!--    Sink or source of humidity due to canopy elements
1843       IF ( plant_canopy ) CALL pcm_tendency( 7 )
1844
1845!
1846!--    Large scale advection. Not implemented for scalars so far.
1847!        IF ( large_scale_forcing )  THEN
1848!           CALL ls_advec( simulated_time, 'q' )
1849!        ENDIF
1850
1851!
1852!--    Nudging. Not implemented for scalars so far.
1853!        IF ( nudging )  CALL nudge( simulated_time, 'q' )
1854
1855!
1856!--    If required compute influence of large-scale subsidence/ascent.
1857!--    Not implemented for scalars so far.
1858       IF ( large_scale_subsidence  .AND.                                      &
1859            .NOT. use_subsidence_tendencies  .AND.                             &
1860            .NOT. large_scale_forcing )  THEN
1861         CALL subsidence( tend, s, s_init, 3 )
1862       ENDIF
1863
1864       CALL user_actions( 's-tendency' )
1865
1866!
1867!--    Prognostic equation for total water content
1868       DO  i = nxl, nxr
1869          DO  j = nys, nyn
1870             DO  k = nzb_s_inner(j,i)+1, nzt
1871                s_p(k,j,i) = s(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +          &
1872                                                  tsc(3) * ts_m(k,j,i) )       &
1873                                      - tsc(5) * rdf_sc(k) *                   &
1874                                        ( s(k,j,i) - s_init(k) )
1875                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
1876             ENDDO
1877          ENDDO
1878       ENDDO
1879
1880!
1881!--    Calculate tendencies for the next Runge-Kutta step
1882       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1883          IF ( intermediate_timestep_count == 1 )  THEN
1884             DO  i = nxl, nxr
1885                DO  j = nys, nyn
1886                   DO  k = nzb_s_inner(j,i)+1, nzt
1887                      ts_m(k,j,i) = tend(k,j,i)
1888                   ENDDO
1889                ENDDO
1890             ENDDO
1891          ELSEIF ( intermediate_timestep_count < &
1892                   intermediate_timestep_count_max )  THEN
1893             DO  i = nxl, nxr
1894                DO  j = nys, nyn
1895                   DO  k = nzb_s_inner(j,i)+1, nzt
1896                      ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * ts_m(k,j,i)
1897                   ENDDO
1898                ENDDO
1899             ENDDO
1900          ENDIF
1901       ENDIF
1902
1903       CALL cpu_log( log_point(66), 's-equation', 'stop' )
1904
1905    ENDIF
1906!
1907!-- If required, compute prognostic equation for turbulent kinetic
1908!-- energy (TKE)
1909    IF ( .NOT. constant_diffusion )  THEN
1910
1911       CALL cpu_log( log_point(16), 'tke-equation', 'start' )
1912
1913       sbt = tsc(2)
1914       IF ( .NOT. use_upstream_for_tke )  THEN
1915          IF ( scalar_advec == 'bc-scheme' )  THEN
1916
1917             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1918!
1919!--             Bott-Chlond scheme always uses Euler time step. Thus:
1920                sbt = 1.0_wp
1921             ENDIF
1922             tend = 0.0_wp
1923             CALL advec_s_bc( e, 'e' )
1924
1925          ENDIF
1926       ENDIF
1927
1928!
1929!--    TKE-tendency terms with no communication
1930       IF ( scalar_advec /= 'bc-scheme'  .OR.  use_upstream_for_tke )  THEN
1931          IF ( use_upstream_for_tke )  THEN
1932             tend = 0.0_wp
1933             CALL advec_s_up( e )
1934          ELSE
1935             tend = 0.0_wp
1936             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1937                IF ( ws_scheme_sca )  THEN
1938                   CALL advec_s_ws( e, 'e' )
1939                ELSE
1940                   CALL advec_s_pw( e )
1941                ENDIF
1942             ELSE
1943                CALL advec_s_up( e )
1944             ENDIF
1945          ENDIF
1946       ENDIF
1947
1948       IF ( .NOT. humidity )  THEN
1949          IF ( ocean )  THEN
1950             CALL diffusion_e( prho, prho_reference )
1951          ELSE
1952             CALL diffusion_e( pt, pt_reference )
1953          ENDIF
1954       ELSE
1955          CALL diffusion_e( vpt, pt_reference )
1956       ENDIF
1957
1958       CALL production_e
1959
1960!
1961!--    Additional sink term for flows through plant canopies
1962       IF ( plant_canopy )  CALL pcm_tendency( 6 )
1963       CALL user_actions( 'e-tendency' )
1964
1965!
1966!--    Prognostic equation for TKE.
1967!--    Eliminate negative TKE values, which can occur due to numerical
1968!--    reasons in the course of the integration. In such cases the old TKE
1969!--    value is reduced by 90%.
1970       DO  i = nxl, nxr
1971          DO  j = nys, nyn
1972             DO  k = nzb_s_inner(j,i)+1, nzt
1973                e_p(k,j,i) = e(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +          &
1974                                                  tsc(3) * te_m(k,j,i) )
1975                IF ( e_p(k,j,i) < 0.0_wp )  e_p(k,j,i) = 0.1_wp * e(k,j,i)
1976             ENDDO
1977          ENDDO
1978       ENDDO
1979
1980!
1981!--    Calculate tendencies for the next Runge-Kutta step
1982       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1983          IF ( intermediate_timestep_count == 1 )  THEN
1984             DO  i = nxl, nxr
1985                DO  j = nys, nyn
1986                   DO  k = nzb_s_inner(j,i)+1, nzt
1987                      te_m(k,j,i) = tend(k,j,i)
1988                   ENDDO
1989                ENDDO
1990             ENDDO
1991          ELSEIF ( intermediate_timestep_count < &
1992                   intermediate_timestep_count_max )  THEN
1993             DO  i = nxl, nxr
1994                DO  j = nys, nyn
1995                   DO  k = nzb_s_inner(j,i)+1, nzt
1996                      te_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * te_m(k,j,i)
1997                   ENDDO
1998                ENDDO
1999             ENDDO
2000          ENDIF
2001       ENDIF
2002
2003       CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
2004
2005    ENDIF
2006
2007 END SUBROUTINE prognostic_equations_vector
2008
2009
2010 END MODULE prognostic_equations_mod
Note: See TracBrowser for help on using the repository browser.