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

Last change on this file since 2809 was 2766, checked in by kanani, 7 years ago

Removal of chem directive, plus minor changes

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