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

Last change on this file since 2118 was 2118, checked in by raasch, 7 years ago

all OpenACC directives and related parts removed from the code

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