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

Last change on this file since 3021 was 3021, checked in by maronga, 3 years ago

bugfixes for nested runs

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