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

Last change on this file since 3274 was 3274, checked in by knoop, 3 years ago

Modularization of all bulk cloud physics code components

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