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

Last change on this file since 3294 was 3294, checked in by raasch, 3 years ago

modularization of the ocean code

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