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

Last change on this file since 3183 was 3183, checked in by suehring, 3 years ago

last commit documented

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