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

Last change on this file since 2233 was 2233, checked in by suehring, 7 years ago

last commit documented

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