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

Last change on this file since 2852 was 2815, checked in by kanani, 6 years ago

Enable restarts and vector optimization for chemistry model

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