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

Last change on this file since 2746 was 2746, checked in by suehring, 4 years ago

Read information from statitic driver for resolved vegetation independently from land- or urban-surface model

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