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

Last change on this file since 4109 was 4109, checked in by suehring, 2 years ago

Control discretization of advection term: separate initialization of WS advection flags for momentum and scalars. In this context, resort the bits and do some minor formatting; Make initialization of scalar-advection flags more flexible, i.e. introduce an arguemnt list to indicate non-cyclic boundaries (required for decycled scalars such as chemical species or aerosols); Introduce extended 'degradation zones', where horizontal advection of passive scalars is discretized by first-order scheme at all grid points that in the vicinity of buildings (<= 3 grid points). Even though no building is within the numerical stencil, first-order scheme is used. At fourth and fifth grid point the order of the horizontal advection scheme is successively upgraded. These extended degradation zones are used to avoid stationary numerical oscillations, which are responsible for high concentration maxima that may appear under shear-free stable conditions. Therefore, an additional 3D interger array used to store control flags is introduced; Change interface for scalar advection routine; Bugfix, avoid uninitialized value sk_num in vector version of WS scalar advection; Chemistry: Decycling boundary conditions are only set at the ghost points not on the prognostic grid points; Land-surface model: Relax checks for non-consistent initialization in case static or dynamic input is provided. For example, soil_temperature or deep_soil_temperature is not mandatory any more if dynamic input is available. Also, improper settings of x_type in namelist are only checked if no static file is available.

  • Property svn:keywords set to Id
File size: 67.7 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-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22! pass integer flag array to WS scalar advection routine which is now necessary
23! as the flags may differ for scalars, e.g. pt can be cyclic while chemical
24! species may be non-cyclic. Further, pass boundary flags.
25!
26! Former revisions:
27! -----------------
28! $Id: prognostic_equations.f90 4109 2019-07-22 17:00:34Z suehring $
29! Application of monotonic flux limiter for the vertical scalar advection
30! up to the topography top (only for the cache-optimized version at the
31! moment). Please note, at the moment the limiter is only applied for passive
32! scalars.
33!
34! 4048 2019-06-21 21:00:21Z knoop
35! Moved tcm_prognostic_equations to module_interface
36!
37! 3987 2019-05-22 09:52:13Z kanani
38! Introduce alternative switch for debug output during timestepping
39!
40! 3956 2019-05-07 12:32:52Z monakurppa
41! Removed salsa calls.
42!
43! 3931 2019-04-24 16:34:28Z schwenkel
44! Correct/complete module_interface introduction for chemistry model
45!
46! 3899 2019-04-16 14:05:27Z monakurppa
47! Corrections in the OpenMP version of salsa
48!
49! 3887 2019 -04-12 08:47:41Z schwenkel
50! Implicit Bugfix for chemistry model, loop for non_transport_physics over
51! ghost points is avoided. Instead introducing module_interface_exchange_horiz.
52!
53! 3885 2019-04-11 11:29:34Z kanani
54! Changes related to global restructuring of location messages and introduction
55! of additional debug messages
56!
57! 3881 2019-04-10 09:31:22Z suehring
58! Bugfix in OpenMP directive
59!
60! 3880 2019-04-08 21:43:02Z knoop
61! Moved wtm_tendencies to module_interface_actions
62!
63! 3874 2019-04-08 16:53:48Z knoop
64! Added non_transport_physics module interfaces and moved bcm code into it
65!
66! 3872 2019-04-08 15:03:06Z knoop
67! Moving prognostic equations of bcm into bulk_cloud_model_mod
68!
69! 3864 2019-04-05 09:01:56Z monakurppa
70! Modifications made for salsa:
71! - salsa_prognostic_equations moved to salsa_mod (and the call to
72!   module_interface_mod)
73! - Renamed nbins --> nbins_aerosol, ncc_tot --> ncomponents_mass and
74!   ngast --> ngases_salsa and loop indices b, c and sg to ib, ic and ig
75!
76! 3840 2019-03-29 10:35:52Z knoop
77! added USE chem_gasphase_mod for nspec, nspec and spc_names
78!
79! 3820 2019-03-27 11:53:41Z forkel
80! renamed do_depo to deposition_dry (ecc)
81!
82! 3797 2019-03-15 11:15:38Z forkel
83! Call chem_integegrate in OpenMP loop   (ketelsen)
84!
85!
86! 3771 2019-02-28 12:19:33Z raasch
87! preprocessor directivs fro rrtmg added
88!
89! 3761 2019-02-25 15:31:42Z raasch
90! unused variable removed
91!
92! 3719 2019-02-06 13:10:18Z kanani
93! Cleaned up chemistry cpu measurements
94!
95! 3684 2019-01-20 20:20:58Z knoop
96! OpenACC port for SPEC
97!
98! 3589 2018-11-30 15:09:51Z suehring
99! Move the control parameter "salsa" from salsa_mod to control_parameters
100! (M. Kurppa)
101!
102! 3582 2018-11-29 19:16:36Z suehring
103! Implementation of a new aerosol module salsa.
104! Remove cpu-logs from i,j loop in cache version.
105!
106! 3458 2018-10-30 14:51:23Z kanani
107! remove duplicate USE chem_modules
108! from chemistry branch r3443, banzhafs, basit:
109! chem_depo call introduced
110! code added for decycling chemistry
111!
112! 3386 2018-10-19 16:28:22Z gronemeier
113! Renamed tcm_prognostic to tcm_prognostic_equations
114!
115! 3355 2018-10-16 14:03:34Z knoop
116! (from branch resler)
117! Fix for chemistry call
118!
119! 3302 2018-10-03 02:39:40Z raasch
120! Stokes drift + wave breaking term added
121!
122! 3298 2018-10-02 12:21:11Z kanani
123! Code added for decycling chemistry (basit)
124!
125! 3294 2018-10-01 02:37:10Z raasch
126! changes concerning modularization of ocean option
127!
128! 3274 2018-09-24 15:42:55Z knoop
129! Modularization of all bulk cloud physics code components
130!
131! 3241 2018-09-12 15:02:00Z raasch
132! omp_get_thread_num now declared in omp directive
133!
134! 3183 2018-07-27 14:25:55Z suehring
135! Remove unused variables from USE statements
136!
137! 3182 2018-07-27 13:36:03Z suehring
138! Revise recent bugfix for nesting
139!
140! 3021 2018-05-16 08:14:20Z maronga
141! Bugfix in IF clause for nesting
142!
143! 3014 2018-05-09 08:42:38Z maronga
144! Fixed a bug in the IF condition to call pcm_tendency in case of
145! potential temperature
146!
147! 2815 2018-02-19 11:29:57Z kanani
148! Rename chem_tendency to chem_prognostic_equations,
149! implement vector version for air chemistry
150!
151! 2766 2018-01-22 17:17:47Z kanani
152! Removed preprocessor directive __chem
153!
154! 2746 2018-01-15 12:06:04Z suehring
155! Move flag plant canopy to modules
156!
157! 2719 2018-01-02 09:02:06Z maronga
158! Bugfix for last change.
159!
160! 2718 2018-01-02 08:49:38Z maronga
161! Corrected "Former revisions" section
162!
163! 2696 2017-12-14 17:12:51Z kanani
164! - Change in file header (GPL part)
165! - Moved TKE equation to tcm_prognostic (TG)
166! - Added switch for chemical reactions (RF, FK)
167! - Implementation of chemistry module (RF, BK, FK)
168!
169! 2563 2017-10-19 15:36:10Z Giersch
170! Variable wind_turbine moved to module control_parameters
171!
172! 2320 2017-07-21 12:47:43Z suehring
173! Modularize large-scale forcing and nudging
174!
175! 2292 2017-06-20 09:51:42Z schwenkel
176! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
177! includes two more prognostic equations for cloud drop concentration (nc) 
178! and cloud water content (qc).
179!
180! 2261 2017-06-08 14:25:57Z raasch
181! bugfix for r2232: openmp directives removed
182!
183! 2233 2017-05-30 18:08:54Z suehring
184!
185! 2232 2017-05-30 17:47:52Z suehring
186! Adjutst to new surface-type structure. Remove call for usm_wall_heat_flux,
187! which is realized directly in diffusion_s now.
188!
189! 2192 2017-03-22 04:14:10Z raasch
190! Bugfix for misplaced and missing openMP directives from r2155
191!
192! 2155 2017-02-21 09:57:40Z hoffmann
193! Bugfix in the calculation of microphysical quantities on ghost points.
194!
195! 2118 2017-01-17 16:38:49Z raasch
196! OpenACC version of subroutine removed
197!
198! 2031 2016-10-21 15:11:58Z knoop
199! renamed variable rho to rho_ocean
200!
201! 2011 2016-09-19 17:29:57Z kanani
202! Flag urban_surface is now defined in module control_parameters.
203!
204! 2007 2016-08-24 15:47:17Z kanani
205! Added pt tendency calculation based on energy balance at urban surfaces
206! (new urban surface model)
207!
208! 2000 2016-08-20 18:09:15Z knoop
209! Forced header and separation lines into 80 columns
210!
211! 1976 2016-07-27 13:28:04Z maronga
212! Simplied calls to radiation model
213!
214! 1960 2016-07-12 16:34:24Z suehring
215! Separate humidity and passive scalar
216!
217! 1914 2016-05-26 14:44:07Z witha
218! Added calls for wind turbine model
219!
220! 1873 2016-04-18 14:50:06Z maronga
221! Module renamed (removed _mod)
222!
223! 1850 2016-04-08 13:29:27Z maronga
224! Module renamed
225!
226! 1826 2016-04-07 12:01:39Z maronga
227! Renamed canopy model calls.
228!
229! 1822 2016-04-07 07:49:42Z hoffmann
230! Kessler microphysics scheme moved to microphysics.
231!
232! 1757 2016-02-22 15:49:32Z maronga
233!
234! 1691 2015-10-26 16:17:44Z maronga
235! Added optional model spin-up without radiation / land surface model calls.
236! Formatting corrections.
237!
238! 1682 2015-10-07 23:56:08Z knoop
239! Code annotations made doxygen readable
240!
241! 1585 2015-04-30 07:05:52Z maronga
242! Added call for temperature tendency calculation due to radiative flux divergence
243!
244! 1517 2015-01-07 19:12:25Z hoffmann
245! advec_s_bc_mod addded, since advec_s_bc is now a module
246!
247! 1484 2014-10-21 10:53:05Z kanani
248! Changes due to new module structure of the plant canopy model:
249! parameters cthf and plant_canopy moved to module plant_canopy_model_mod.
250! Removed double-listing of use_upstream_for_tke in ONLY-list of module
251! control_parameters
252!
253! 1409 2014-05-23 12:11:32Z suehring
254! Bugfix: i_omp_start changed for advec_u_ws at left inflow and outflow boundary.
255! This ensures that left-hand side fluxes are also calculated for nxl in that
256! case, even though the solution at nxl is overwritten in boundary_conds()
257!
258! 1398 2014-05-07 11:15:00Z heinze
259! Rayleigh-damping for horizontal velocity components changed: instead of damping
260! against ug and vg, damping against u_init and v_init is used to allow for a
261! homogenized treatment in case of nudging
262!
263! 1380 2014-04-28 12:40:45Z heinze
264! Change order of calls for scalar prognostic quantities:
265! ls_advec -> nudging -> subsidence since initial profiles
266!
267! 1374 2014-04-25 12:55:07Z raasch
268! missing variables added to ONLY lists
269!
270! 1365 2014-04-22 15:03:56Z boeske
271! Calls of ls_advec for large scale advection added,
272! subroutine subsidence is only called if use_subsidence_tendencies = .F.,
273! new argument ls_index added to the calls of subsidence
274! +ls_index
275!
276! 1361 2014-04-16 15:17:48Z hoffmann
277! Two-moment microphysics moved to the start of prognostic equations. This makes
278! the 3d arrays for tend_q, tend_qr, tend_pt and tend_pt redundant.
279! Additionally, it is allowed to call the microphysics just once during the time
280! step (not at each sub-time step).
281!
282! Two-moment cloud physics added for vector and accelerator optimization.
283!
284! 1353 2014-04-08 15:21:23Z heinze
285! REAL constants provided with KIND-attribute
286!
287! 1337 2014-03-25 15:11:48Z heinze
288! Bugfix: REAL constants provided with KIND-attribute
289!
290! 1332 2014-03-25 11:59:43Z suehring
291! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke
292!
293! 1330 2014-03-24 17:29:32Z suehring
294! In case of SGS-particle velocity advection of TKE is also allowed with
295! dissipative 5th-order scheme.
296!
297! 1320 2014-03-20 08:40:49Z raasch
298! ONLY-attribute added to USE-statements,
299! kind-parameters added to all INTEGER and REAL declaration statements,
300! kinds are defined in new module kinds,
301! old module precision_kind is removed,
302! revision history before 2012 removed,
303! comment fields (!:) to be used for variable explanations added to
304! all variable declaration statements
305!
306! 1318 2014-03-17 13:35:16Z raasch
307! module interfaces removed
308!
309! 1257 2013-11-08 15:18:40Z raasch
310! openacc loop vector clauses removed, independent clauses added
311!
312! 1246 2013-11-01 08:59:45Z heinze
313! enable nudging also for accelerator version
314!
315! 1241 2013-10-30 11:36:58Z heinze
316! usage of nudging enabled (so far not implemented for accelerator version)
317!
318! 1179 2013-06-14 05:57:58Z raasch
319! two arguments removed from routine buoyancy, ref_state updated on device
320!
321! 1128 2013-04-12 06:19:32Z raasch
322! those parts requiring global communication moved to time_integration,
323! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
324! j_north
325!
326! 1115 2013-03-26 18:16:16Z hoffmann
327! optimized cloud physics: calculation of microphysical tendencies transfered
328! to microphysics.f90; qr and nr are only calculated if precipitation is required
329!
330! 1111 2013-03-08 23:54:10Z raasch
331! update directives for prognostic quantities removed
332!
333! 1106 2013-03-04 05:31:38Z raasch
334! small changes in code formatting
335!
336! 1092 2013-02-02 11:24:22Z raasch
337! unused variables removed
338!
339! 1053 2012-11-13 17:11:03Z hoffmann
340! implementation of two new prognostic equations for rain drop concentration (nr)
341! and rain water content (qr)
342!
343! currently, only available for cache loop optimization
344!
345! 1036 2012-10-22 13:43:42Z raasch
346! code put under GPL (PALM 3.9)
347!
348! 1019 2012-09-28 06:46:45Z raasch
349! non-optimized version of prognostic_equations removed
350!
351! 1015 2012-09-27 09:23:24Z raasch
352! new branch prognostic_equations_acc
353! OpenACC statements added + code changes required for GPU optimization
354!
355! 1001 2012-09-13 14:08:46Z raasch
356! all actions concerning leapfrog- and upstream-spline-scheme removed
357!
358! 978 2012-08-09 08:28:32Z fricke
359! km_damp_x and km_damp_y removed in calls of diffusion_u and diffusion_v
360! add ptdf_x, ptdf_y for damping the potential temperature at the inflow
361! boundary in case of non-cyclic lateral boundaries
362! Bugfix: first thread index changes for WS-scheme at the inflow
363!
364! 940 2012-07-09 14:31:00Z raasch
365! temperature equation can be switched off
366!
367! Revision 1.1  2000/04/13 14:56:27  schroeter
368! Initial revision
369!
370!
371! Description:
372! ------------
373!> Solving the prognostic equations.
374!------------------------------------------------------------------------------!
375 MODULE prognostic_equations_mod
376
377    USE advec_s_bc_mod,                                                        &
378        ONLY:  advec_s_bc
379
380    USE advec_s_pw_mod,                                                        &
381        ONLY:  advec_s_pw
382
383    USE advec_s_up_mod,                                                        &
384        ONLY:  advec_s_up
385
386    USE advec_u_pw_mod,                                                        &
387        ONLY:  advec_u_pw
388
389    USE advec_u_up_mod,                                                        &
390        ONLY:  advec_u_up
391
392    USE advec_v_pw_mod,                                                        &
393        ONLY:  advec_v_pw
394
395    USE advec_v_up_mod,                                                        &
396        ONLY:  advec_v_up
397
398    USE advec_w_pw_mod,                                                        &
399        ONLY:  advec_w_pw
400
401    USE advec_w_up_mod,                                                        &
402        ONLY:  advec_w_up
403
404    USE advec_ws,                                                              &
405        ONLY:  advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws
406
407    USE arrays_3d,                                                             &
408        ONLY:  diss_l_e, diss_l_pt, diss_l_q,                                  &
409               diss_l_s, diss_l_sa, diss_s_e,                                  &
410               diss_s_pt, diss_s_q, diss_s_s, diss_s_sa,                       &
411               e, e_p, flux_s_e, flux_s_pt, flux_s_q,                          &
412               flux_s_s, flux_s_sa, flux_l_e,                                  &
413               flux_l_pt, flux_l_q, flux_l_s,                                  &
414               flux_l_sa, pt, ptdf_x, ptdf_y, pt_init,                         &
415               pt_p, prho, q, q_init, q_p, rdf, rdf_sc,                        &
416               ref_state, rho_ocean, s, s_init, s_p, tend, te_m,               &
417               tpt_m, tq_m, ts_m, tu_m, tv_m, tw_m, u,                         &
418               ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
419
420    USE buoyancy_mod,                                                          &
421        ONLY:  buoyancy
422
423    USE control_parameters,                                                    &
424        ONLY:  bc_dirichlet_l,                                                 &
425               bc_dirichlet_n,                                                 &
426               bc_dirichlet_r,                                                 &
427               bc_dirichlet_s,                                                 &
428               bc_radiation_l,                                                 &
429               bc_radiation_n,                                                 &
430               bc_radiation_r,                                                 &
431               bc_radiation_s,                                                 &
432               constant_diffusion,                                             &
433               debug_output_timestep,                                          &
434               dp_external, dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d,    &
435               humidity, intermediate_timestep_count,                          &
436               intermediate_timestep_count_max, large_scale_forcing,           &
437               large_scale_subsidence,                                         &
438               monotonic_limiter_z,                                            &
439               neutral, nudging,                                               &
440               ocean_mode, passive_scalar, plant_canopy, pt_reference,         &
441               scalar_advec, scalar_advec, simulated_time, sloping_surface,    &
442               timestep_scheme, tsc, use_subsidence_tendencies,                &
443               use_upstream_for_tke, wind_turbine, ws_scheme_mom,              & 
444               ws_scheme_sca, urban_surface, land_surface,                     &
445               time_since_reference_point, salsa
446
447    USE coriolis_mod,                                                          &
448        ONLY:  coriolis
449
450    USE cpulog,                                                                &
451        ONLY:  cpu_log, log_point, log_point_s
452
453    USE diffusion_s_mod,                                                       &
454        ONLY:  diffusion_s
455
456    USE diffusion_u_mod,                                                       &
457        ONLY:  diffusion_u
458
459    USE diffusion_v_mod,                                                       &
460        ONLY:  diffusion_v
461
462    USE diffusion_w_mod,                                                       &
463        ONLY:  diffusion_w
464
465    USE indices,                                                               &
466        ONLY:  advc_flags_s,                                                   &
467               nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
468               nzb, nzt, wall_flags_0
469
470    USE kinds
471
472    USE lsf_nudging_mod,                                                       &
473        ONLY:  ls_advec, nudge
474
475    USE module_interface,                                                      &
476        ONLY:  module_interface_actions, &
477               module_interface_non_advective_processes, &
478               module_interface_exchange_horiz, &
479               module_interface_prognostic_equations
480
481    USE ocean_mod,                                                             &
482        ONLY:  stokes_drift_terms, stokes_force,   &
483               wave_breaking, wave_breaking_term
484
485    USE plant_canopy_model_mod,                                                &
486        ONLY:  cthf, pcm_tendency
487
488#if defined( __rrtmg )
489    USE radiation_model_mod,                                                   &
490        ONLY:  radiation, radiation_tendency,                                  &
491               skip_time_do_radiation
492#endif
493
494    USE statistics,                                                            &
495        ONLY:  hom
496
497    USE subsidence_mod,                                                        &
498        ONLY:  subsidence
499
500    USE surface_mod,                                                           &
501        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
502                surf_usm_v
503
504    IMPLICIT NONE
505
506    PRIVATE
507    PUBLIC prognostic_equations_cache, prognostic_equations_vector
508
509    INTERFACE prognostic_equations_cache
510       MODULE PROCEDURE prognostic_equations_cache
511    END INTERFACE prognostic_equations_cache
512
513    INTERFACE prognostic_equations_vector
514       MODULE PROCEDURE prognostic_equations_vector
515    END INTERFACE prognostic_equations_vector
516
517
518 CONTAINS
519
520
521!------------------------------------------------------------------------------!
522! Description:
523! ------------
524!> Version with one optimized loop over all equations. It is only allowed to
525!> be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
526!>
527!> Here the calls of most subroutines are embedded in two DO loops over i and j,
528!> so communication between CPUs is not allowed (does not make sense) within
529!> these loops.
530!>
531!> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
532!------------------------------------------------------------------------------!
533
534 SUBROUTINE prognostic_equations_cache
535
536
537    INTEGER(iwp) ::  i                   !<
538    INTEGER(iwp) ::  i_omp_start         !<
539    INTEGER(iwp) ::  j                   !<
540    INTEGER(iwp) ::  k                   !<
541!$  INTEGER(iwp) ::  omp_get_thread_num  !<
542    INTEGER(iwp) ::  tn = 0              !<
543
544    LOGICAL      ::  loop_start          !<
545
546
547    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_cache', 'start' )
548!
549!-- Time measurement can only be performed for the whole set of equations
550    CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
551
552    !$OMP PARALLEL PRIVATE (i,j)
553    !$OMP DO
554    DO  i = nxl, nxr
555       DO  j = nys, nyn
556!
557!--       Calculate non advective processes for all other modules
558          CALL module_interface_non_advective_processes( i, j )
559       ENDDO
560    ENDDO
561!
562!-- Module Inferface for exchange horiz after non_advective_processes but before
563!-- advection. Therefore, non_advective_processes must not run for ghost points.
564    !$OMP END PARALLEL
565    CALL module_interface_exchange_horiz()
566!
567!-- Loop over all prognostic equations
568    !$OMP PARALLEL PRIVATE (i,i_omp_start,j,k,loop_start,tn)
569
570    !$  tn = omp_get_thread_num()
571    loop_start = .TRUE.
572
573    !$OMP DO
574    DO  i = nxl, nxr
575
576!
577!--    Store the first loop index. It differs for each thread and is required
578!--    later in advec_ws
579       IF ( loop_start )  THEN
580          loop_start  = .FALSE.
581          i_omp_start = i
582       ENDIF
583
584       DO  j = nys, nyn
585!
586!--       Tendency terms for u-velocity component. Please note, in case of
587!--       non-cyclic boundary conditions the grid point i=0 is excluded from
588!--       the prognostic equations for the u-component.
589          IF ( i >= nxlu )  THEN
590
591             tend(:,j,i) = 0.0_wp
592             IF ( timestep_scheme(1:5) == 'runge' )  THEN
593                IF ( ws_scheme_mom )  THEN
594                   CALL advec_u_ws( i, j, i_omp_start, tn )
595                ELSE
596                   CALL advec_u_pw( i, j )
597                ENDIF
598             ELSE
599                CALL advec_u_up( i, j )
600             ENDIF
601             CALL diffusion_u( i, j )
602             CALL coriolis( i, j, 1 )
603             IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
604                CALL buoyancy( i, j, pt, 1 )
605             ENDIF
606
607!
608!--          Drag by plant canopy
609             IF ( plant_canopy )  CALL pcm_tendency( i, j, 1 )
610
611!
612!--          External pressure gradient
613             IF ( dp_external )  THEN
614                DO  k = dp_level_ind_b+1, nzt
615                   tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
616                ENDDO
617             ENDIF
618
619!
620!--          Nudging
621             IF ( nudging )  CALL nudge( i, j, simulated_time, 'u' )
622
623!
624!--          Effect of Stokes drift (in ocean mode only)
625             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 1 )
626
627             CALL module_interface_actions( i, j, 'u-tendency' )
628!
629!--          Prognostic equation for u-velocity component
630             DO  k = nzb+1, nzt
631
632                u_p(k,j,i) = u(k,j,i) + ( dt_3d *                               &
633                                            ( tsc(2) * tend(k,j,i) +            &
634                                              tsc(3) * tu_m(k,j,i) )            &
635                                            - tsc(5) * rdf(k)                   &
636                                                     * ( u(k,j,i) - u_init(k) ) &
637                                        ) * MERGE( 1.0_wp, 0.0_wp,              &
638                                                   BTEST( wall_flags_0(k,j,i), 1 )&
639                                                 )
640             ENDDO
641
642!
643!--          Add turbulence generated by wave breaking (in ocean mode only)
644             IF ( wave_breaking  .AND.                                         &
645               intermediate_timestep_count == intermediate_timestep_count_max )&
646             THEN
647                CALL wave_breaking_term( i, j, 1 )
648             ENDIF
649
650!
651!--          Calculate tendencies for the next Runge-Kutta step
652             IF ( timestep_scheme(1:5) == 'runge' )  THEN
653                IF ( intermediate_timestep_count == 1 )  THEN
654                   DO  k = nzb+1, nzt
655                      tu_m(k,j,i) = tend(k,j,i)
656                   ENDDO
657                ELSEIF ( intermediate_timestep_count < &
658                         intermediate_timestep_count_max )  THEN
659                   DO  k = nzb+1, nzt
660                      tu_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                &
661                                     + 5.3125_wp * tu_m(k,j,i)
662                   ENDDO
663                ENDIF
664             ENDIF
665
666          ENDIF
667!
668!--       Tendency terms for v-velocity component. Please note, in case of
669!--       non-cyclic boundary conditions the grid point j=0 is excluded from
670!--       the prognostic equations for the v-component. !--       
671          IF ( j >= nysv )  THEN
672
673             tend(:,j,i) = 0.0_wp
674             IF ( timestep_scheme(1:5) == 'runge' )  THEN
675                IF ( ws_scheme_mom )  THEN
676                    CALL advec_v_ws( i, j, i_omp_start, tn )
677                ELSE
678                    CALL advec_v_pw( i, j )
679                ENDIF
680             ELSE
681                CALL advec_v_up( i, j )
682             ENDIF
683             CALL diffusion_v( i, j )
684             CALL coriolis( i, j, 2 )
685
686!
687!--          Drag by plant canopy
688             IF ( plant_canopy )  CALL pcm_tendency( i, j, 2 )
689
690!
691!--          External pressure gradient
692             IF ( dp_external )  THEN
693                DO  k = dp_level_ind_b+1, nzt
694                   tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
695                ENDDO
696             ENDIF
697
698!
699!--          Nudging
700             IF ( nudging )  CALL nudge( i, j, simulated_time, 'v' )
701
702!
703!--          Effect of Stokes drift (in ocean mode only)
704             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 2 )
705
706             CALL module_interface_actions( i, j, 'v-tendency' )
707!
708!--          Prognostic equation for v-velocity component
709             DO  k = nzb+1, nzt
710                v_p(k,j,i) = v(k,j,i) + ( dt_3d *                              &
711                                            ( tsc(2) * tend(k,j,i) +           &
712                                              tsc(3) * tv_m(k,j,i) )           &
713                                            - tsc(5) * rdf(k)                  &
714                                                     * ( v(k,j,i) - v_init(k) )&
715                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
716                                                   BTEST( wall_flags_0(k,j,i), 2 )&
717                                                 )
718             ENDDO
719
720!
721!--          Add turbulence generated by wave breaking (in ocean mode only)
722             IF ( wave_breaking  .AND.                                         &
723               intermediate_timestep_count == intermediate_timestep_count_max )&
724             THEN
725                CALL wave_breaking_term( i, j, 2 )
726             ENDIF
727
728!
729!--          Calculate tendencies for the next Runge-Kutta step
730             IF ( timestep_scheme(1:5) == 'runge' )  THEN
731                IF ( intermediate_timestep_count == 1 )  THEN
732                   DO  k = nzb+1, nzt
733                      tv_m(k,j,i) = tend(k,j,i)
734                   ENDDO
735                ELSEIF ( intermediate_timestep_count < &
736                         intermediate_timestep_count_max )  THEN
737                   DO  k = nzb+1, nzt
738                      tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
739                                     + 5.3125_wp * tv_m(k,j,i)
740                   ENDDO
741                ENDIF
742             ENDIF
743
744          ENDIF
745
746!
747!--       Tendency terms for w-velocity component
748          tend(:,j,i) = 0.0_wp
749          IF ( timestep_scheme(1:5) == 'runge' )  THEN
750             IF ( ws_scheme_mom )  THEN
751                CALL advec_w_ws( i, j, i_omp_start, tn )
752             ELSE
753                CALL advec_w_pw( i, j )
754             END IF
755          ELSE
756             CALL advec_w_up( i, j )
757          ENDIF
758          CALL diffusion_w( i, j )
759          CALL coriolis( i, j, 3 )
760
761          IF ( .NOT. neutral )  THEN
762             IF ( ocean_mode )  THEN
763                CALL buoyancy( i, j, rho_ocean, 3 )
764             ELSE
765                IF ( .NOT. humidity )  THEN
766                   CALL buoyancy( i, j, pt, 3 )
767                ELSE
768                   CALL buoyancy( i, j, vpt, 3 )
769                ENDIF
770             ENDIF
771          ENDIF
772
773!
774!--       Drag by plant canopy
775          IF ( plant_canopy )  CALL pcm_tendency( i, j, 3 )
776
777!
778!--       Effect of Stokes drift (in ocean mode only)
779          IF ( stokes_force )  CALL stokes_drift_terms( i, j, 3 )
780
781          CALL module_interface_actions( i, j, 'w-tendency' )
782!
783!--       Prognostic equation for w-velocity component
784          DO  k = nzb+1, nzt-1
785             w_p(k,j,i) = w(k,j,i) + ( dt_3d *                                 &
786                                             ( tsc(2) * tend(k,j,i) +          &
787                                               tsc(3) * tw_m(k,j,i) )          &
788                                             - tsc(5) * rdf(k) * w(k,j,i)      &
789                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
790                                                BTEST( wall_flags_0(k,j,i), 3 )&
791                                              )
792          ENDDO
793
794!
795!--       Calculate tendencies for the next Runge-Kutta step
796          IF ( timestep_scheme(1:5) == 'runge' )  THEN
797             IF ( intermediate_timestep_count == 1 )  THEN
798                DO  k = nzb+1, nzt-1
799                   tw_m(k,j,i) = tend(k,j,i)
800                ENDDO
801             ELSEIF ( intermediate_timestep_count < &
802                      intermediate_timestep_count_max )  THEN
803                DO  k = nzb+1, nzt-1
804                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
805                                  + 5.3125_wp * tw_m(k,j,i)
806                ENDDO
807             ENDIF
808          ENDIF
809
810!
811!--       If required, compute prognostic equation for potential temperature
812          IF ( .NOT. neutral )  THEN
813!
814!--          Tendency terms for potential temperature
815             tend(:,j,i) = 0.0_wp
816             IF ( timestep_scheme(1:5) == 'runge' )  THEN
817                   IF ( ws_scheme_sca )  THEN
818                      CALL advec_s_ws( advc_flags_s,                           &
819                                       i, j, pt, 'pt', flux_s_pt, diss_s_pt,   &
820                                       flux_l_pt, diss_l_pt, i_omp_start, tn,  &
821                                       bc_dirichlet_l  .OR.  bc_radiation_l,   &
822                                       bc_dirichlet_n  .OR.  bc_radiation_n,   &
823                                       bc_dirichlet_r  .OR.  bc_radiation_r,   &
824                                       bc_dirichlet_s  .OR.  bc_radiation_s )
825                   ELSE
826                      CALL advec_s_pw( i, j, pt )
827                   ENDIF
828             ELSE
829                CALL advec_s_up( i, j, pt )
830             ENDIF
831             CALL diffusion_s( i, j, pt,                                       &
832                               surf_def_h(0)%shf, surf_def_h(1)%shf,           &
833                               surf_def_h(2)%shf,                              &
834                               surf_lsm_h%shf,    surf_usm_h%shf,              &
835                               surf_def_v(0)%shf, surf_def_v(1)%shf,           &
836                               surf_def_v(2)%shf, surf_def_v(3)%shf,           &
837                               surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,           &
838                               surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,           &
839                               surf_usm_v(0)%shf, surf_usm_v(1)%shf,           &
840                               surf_usm_v(2)%shf, surf_usm_v(3)%shf )
841
842!
843!--          Consideration of heat sources within the plant canopy
844             IF ( plant_canopy  .AND.                                          &
845                (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
846                CALL pcm_tendency( i, j, 4 )
847             ENDIF
848
849!
850!--          Large scale advection
851             IF ( large_scale_forcing )  THEN
852                CALL ls_advec( i, j, simulated_time, 'pt' )
853             ENDIF
854
855!
856!--          Nudging
857             IF ( nudging )  CALL nudge( i, j, simulated_time, 'pt' )
858
859!
860!--          If required, compute effect of large-scale subsidence/ascent
861             IF ( large_scale_subsidence  .AND.                                &
862                  .NOT. use_subsidence_tendencies )  THEN
863                CALL subsidence( i, j, tend, pt, pt_init, 2 )
864             ENDIF
865
866#if defined( __rrtmg )
867!
868!--          If required, add tendency due to radiative heating/cooling
869             IF ( radiation  .AND.                                             &
870                  simulated_time > skip_time_do_radiation )  THEN
871                CALL radiation_tendency ( i, j, tend )
872             ENDIF
873#endif
874
875             CALL module_interface_actions( i, j, 'pt-tendency' )
876!
877!--          Prognostic equation for potential temperature
878             DO  k = nzb+1, nzt
879                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d *                            &
880                                                  ( tsc(2) * tend(k,j,i) +     &
881                                                    tsc(3) * tpt_m(k,j,i) )    &
882                                                  - tsc(5)                     &
883                                                  * ( pt(k,j,i) - pt_init(k) ) &
884                                                  * ( rdf_sc(k) + ptdf_x(i)    &
885                                                                + ptdf_y(j) )  &
886                                          )                                    &
887                                       * MERGE( 1.0_wp, 0.0_wp,                &
888                                                BTEST( wall_flags_0(k,j,i), 0 )&
889                                              )
890             ENDDO
891
892!
893!--          Calculate tendencies for the next Runge-Kutta step
894             IF ( timestep_scheme(1:5) == 'runge' )  THEN
895                IF ( intermediate_timestep_count == 1 )  THEN
896                   DO  k = nzb+1, nzt
897                      tpt_m(k,j,i) = tend(k,j,i)
898                   ENDDO
899                ELSEIF ( intermediate_timestep_count < &
900                         intermediate_timestep_count_max )  THEN
901                   DO  k = nzb+1, nzt
902                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
903                                        5.3125_wp * tpt_m(k,j,i)
904                   ENDDO
905                ENDIF
906             ENDIF
907
908          ENDIF
909
910!
911!--       If required, compute prognostic equation for total water content
912          IF ( humidity )  THEN
913
914!
915!--          Tendency-terms for total water content / scalar
916             tend(:,j,i) = 0.0_wp
917             IF ( timestep_scheme(1:5) == 'runge' )                            &
918             THEN
919                IF ( ws_scheme_sca )  THEN
920                   CALL advec_s_ws( advc_flags_s,                              &
921                                    i, j, q, 'q', flux_s_q,                    &
922                                    diss_s_q, flux_l_q, diss_l_q,              &
923                                    i_omp_start, tn,                           &
924                                    bc_dirichlet_l  .OR.  bc_radiation_l,      &
925                                    bc_dirichlet_n  .OR.  bc_radiation_n,      &
926                                    bc_dirichlet_r  .OR.  bc_radiation_r,      &
927                                    bc_dirichlet_s  .OR.  bc_radiation_s )
928                ELSE
929                   CALL advec_s_pw( i, j, q )
930                ENDIF
931             ELSE
932                CALL advec_s_up( i, j, q )
933             ENDIF
934             CALL diffusion_s( i, j, q,                                        &
935                               surf_def_h(0)%qsws, surf_def_h(1)%qsws,         &
936                               surf_def_h(2)%qsws,                             &
937                               surf_lsm_h%qsws,    surf_usm_h%qsws,            &
938                               surf_def_v(0)%qsws, surf_def_v(1)%qsws,         &
939                               surf_def_v(2)%qsws, surf_def_v(3)%qsws,         &
940                               surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,         &
941                               surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,         &
942                               surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,         &
943                               surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
944
945!
946!--          Sink or source of humidity due to canopy elements
947             IF ( plant_canopy )  CALL pcm_tendency( i, j, 5 )
948
949!
950!--          Large scale advection
951             IF ( large_scale_forcing )  THEN
952                CALL ls_advec( i, j, simulated_time, 'q' )
953             ENDIF
954
955!
956!--          Nudging
957             IF ( nudging )  CALL nudge( i, j, simulated_time, 'q' )
958
959!
960!--          If required compute influence of large-scale subsidence/ascent
961             IF ( large_scale_subsidence  .AND.                                &
962                  .NOT. use_subsidence_tendencies )  THEN
963                CALL subsidence( i, j, tend, q, q_init, 3 )
964             ENDIF
965
966             CALL module_interface_actions( i, j, 'q-tendency' )
967
968!
969!--          Prognostic equation for total water content / scalar
970             DO  k = nzb+1, nzt
971                q_p(k,j,i) = q(k,j,i) + ( dt_3d *                              &
972                                                ( tsc(2) * tend(k,j,i) +       &
973                                                  tsc(3) * tq_m(k,j,i) )       &
974                                                - tsc(5) * rdf_sc(k) *         &
975                                                      ( q(k,j,i) - q_init(k) ) &
976                                        )                                      &
977                                       * MERGE( 1.0_wp, 0.0_wp,                &
978                                                BTEST( wall_flags_0(k,j,i), 0 )&
979                                              )               
980                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
981             ENDDO
982
983!
984!--          Calculate tendencies for the next Runge-Kutta step
985             IF ( timestep_scheme(1:5) == 'runge' )  THEN
986                IF ( intermediate_timestep_count == 1 )  THEN
987                   DO  k = nzb+1, nzt
988                      tq_m(k,j,i) = tend(k,j,i)
989                   ENDDO
990                ELSEIF ( intermediate_timestep_count < &
991                         intermediate_timestep_count_max )  THEN
992                   DO  k = nzb+1, nzt
993                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
994                                       5.3125_wp * tq_m(k,j,i)
995                   ENDDO
996                ENDIF
997             ENDIF
998
999          ENDIF
1000
1001!
1002!--       If required, compute prognostic equation for scalar
1003          IF ( passive_scalar )  THEN
1004!
1005!--          Tendency-terms for total water content / scalar
1006             tend(:,j,i) = 0.0_wp
1007             IF ( timestep_scheme(1:5) == 'runge' )                            &
1008             THEN
1009                IF ( ws_scheme_sca )  THEN
1010!
1011!--                For scalar advection apply monotonic flux limiter near
1012!--                topography.
1013                   CALL advec_s_ws( advc_flags_s,                              &
1014                                    i, j, s, 's', flux_s_s,                    &
1015                                    diss_s_s, flux_l_s, diss_l_s, i_omp_start, &
1016                                    tn,                                        &
1017                                    bc_dirichlet_l  .OR.  bc_radiation_l,      &
1018                                    bc_dirichlet_n  .OR.  bc_radiation_n,      &
1019                                    bc_dirichlet_r  .OR.  bc_radiation_r,      &
1020                                    bc_dirichlet_s  .OR.  bc_radiation_s,      &
1021                                    monotonic_limiter_z )
1022                ELSE
1023                   CALL advec_s_pw( i, j, s )
1024                ENDIF
1025             ELSE
1026                CALL advec_s_up( i, j, s )
1027             ENDIF
1028             CALL diffusion_s( i, j, s,                                        &
1029                               surf_def_h(0)%ssws, surf_def_h(1)%ssws,         &
1030                               surf_def_h(2)%ssws,                             &
1031                               surf_lsm_h%ssws,    surf_usm_h%ssws,            &
1032                               surf_def_v(0)%ssws, surf_def_v(1)%ssws,         &
1033                               surf_def_v(2)%ssws, surf_def_v(3)%ssws,         &
1034                               surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,         &
1035                               surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,         &
1036                               surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,         &
1037                               surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1038
1039!
1040!--          Sink or source of scalar concentration due to canopy elements
1041             IF ( plant_canopy )  CALL pcm_tendency( i, j, 7 )
1042
1043!
1044!--          Large scale advection, still need to be extended for scalars
1045!              IF ( large_scale_forcing )  THEN
1046!                 CALL ls_advec( i, j, simulated_time, 's' )
1047!              ENDIF
1048
1049!
1050!--          Nudging, still need to be extended for scalars
1051!              IF ( nudging )  CALL nudge( i, j, simulated_time, 's' )
1052
1053!
1054!--          If required compute influence of large-scale subsidence/ascent.
1055!--          Note, the last argument is of no meaning in this case, as it is
1056!--          only used in conjunction with large_scale_forcing, which is to
1057!--          date not implemented for scalars.
1058             IF ( large_scale_subsidence  .AND.                                &
1059                  .NOT. use_subsidence_tendencies  .AND.                       &
1060                  .NOT. large_scale_forcing )  THEN
1061                CALL subsidence( i, j, tend, s, s_init, 3 )
1062             ENDIF
1063
1064             CALL module_interface_actions( i, j, 's-tendency' )
1065
1066!
1067!--          Prognostic equation for scalar
1068             DO  k = nzb+1, nzt
1069                s_p(k,j,i) = s(k,j,i) + (  dt_3d *                             &
1070                                            ( tsc(2) * tend(k,j,i) +           &
1071                                              tsc(3) * ts_m(k,j,i) )           &
1072                                            - tsc(5) * rdf_sc(k)               &
1073                                                     * ( s(k,j,i) - s_init(k) )&
1074                                        )                                      &
1075                                       * MERGE( 1.0_wp, 0.0_wp,                &
1076                                                BTEST( wall_flags_0(k,j,i), 0 )&
1077                                              )
1078                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
1079             ENDDO
1080
1081!
1082!--          Calculate tendencies for the next Runge-Kutta step
1083             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1084                IF ( intermediate_timestep_count == 1 )  THEN
1085                   DO  k = nzb+1, nzt
1086                      ts_m(k,j,i) = tend(k,j,i)
1087                   ENDDO
1088                ELSEIF ( intermediate_timestep_count < &
1089                         intermediate_timestep_count_max )  THEN
1090                   DO  k = nzb+1, nzt
1091                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
1092                                       5.3125_wp * ts_m(k,j,i)
1093                   ENDDO
1094                ENDIF
1095             ENDIF
1096
1097          ENDIF
1098!
1099!--       Calculate prognostic equations for all other modules
1100          CALL module_interface_prognostic_equations( i, j, i_omp_start, tn )
1101
1102       ENDDO  ! loop over j
1103    ENDDO  ! loop over i
1104    !$OMP END PARALLEL
1105
1106
1107    CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1108
1109    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_cache', 'end' )
1110
1111 END SUBROUTINE prognostic_equations_cache
1112
1113
1114!------------------------------------------------------------------------------!
1115! Description:
1116! ------------
1117!> Version for vector machines
1118!------------------------------------------------------------------------------!
1119
1120 SUBROUTINE prognostic_equations_vector
1121
1122
1123    INTEGER(iwp) ::  i     !<
1124    INTEGER(iwp) ::  j     !<
1125    INTEGER(iwp) ::  k     !<
1126
1127    REAL(wp)     ::  sbt  !<
1128
1129
1130    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_vector', 'start' )
1131!
1132!-- Calculate non advective processes for all other modules
1133    CALL module_interface_non_advective_processes
1134!
1135!-- Module Inferface for exchange horiz after non_advective_processes but before
1136!-- advection. Therefore, non_advective_processes must not run for ghost points.
1137    CALL module_interface_exchange_horiz()
1138!
1139!-- u-velocity component
1140    CALL cpu_log( log_point(5), 'u-equation', 'start' )
1141
1142    !$ACC KERNELS PRESENT(tend)
1143    tend = 0.0_wp
1144    !$ACC END KERNELS
1145    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1146       IF ( ws_scheme_mom )  THEN
1147          CALL advec_u_ws
1148       ELSE
1149          CALL advec_u_pw
1150       ENDIF
1151    ELSE
1152       CALL advec_u_up
1153    ENDIF
1154    CALL diffusion_u
1155    CALL coriolis( 1 )
1156    IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
1157       CALL buoyancy( pt, 1 )
1158    ENDIF
1159
1160!
1161!-- Drag by plant canopy
1162    IF ( plant_canopy )  CALL pcm_tendency( 1 )
1163
1164!
1165!-- External pressure gradient
1166    IF ( dp_external )  THEN
1167       DO  i = nxlu, nxr
1168          DO  j = nys, nyn
1169             DO  k = dp_level_ind_b+1, nzt
1170                tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
1171             ENDDO
1172          ENDDO
1173       ENDDO
1174    ENDIF
1175
1176!
1177!-- Nudging
1178    IF ( nudging )  CALL nudge( simulated_time, 'u' )
1179
1180!
1181!-- Effect of Stokes drift (in ocean mode only)
1182    IF ( stokes_force )  CALL stokes_drift_terms( 1 )
1183
1184    CALL module_interface_actions( 'u-tendency' )
1185
1186!
1187!-- Prognostic equation for u-velocity component
1188    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1189    !$ACC PRESENT(u, tend, tu_m, u_init, rdf, wall_flags_0) &
1190    !$ACC PRESENT(tsc(2:5)) &
1191    !$ACC PRESENT(u_p)
1192    DO  i = nxlu, nxr
1193       DO  j = nys, nyn
1194          DO  k = nzb+1, nzt
1195             u_p(k,j,i) = u(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +          &
1196                                                 tsc(3) * tu_m(k,j,i) )          &
1197                                               - tsc(5) * rdf(k) *               &
1198                                                        ( u(k,j,i) - u_init(k) ) &
1199                                     ) * MERGE( 1.0_wp, 0.0_wp,                  &
1200                                                 BTEST( wall_flags_0(k,j,i), 1 ) &
1201                                              )
1202          ENDDO
1203       ENDDO
1204    ENDDO
1205
1206!
1207!-- Add turbulence generated by wave breaking (in ocean mode only)
1208    IF ( wave_breaking  .AND.                                                  &
1209         intermediate_timestep_count == intermediate_timestep_count_max )      &
1210    THEN
1211       CALL wave_breaking_term( 1 )
1212    ENDIF
1213
1214!
1215!-- Calculate tendencies for the next Runge-Kutta step
1216    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1217       IF ( intermediate_timestep_count == 1 )  THEN
1218          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1219          !$ACC PRESENT(tend, tu_m)
1220          DO  i = nxlu, nxr
1221             DO  j = nys, nyn
1222                DO  k = nzb+1, nzt
1223                   tu_m(k,j,i) = tend(k,j,i)
1224                ENDDO
1225             ENDDO
1226          ENDDO
1227       ELSEIF ( intermediate_timestep_count < &
1228                intermediate_timestep_count_max )  THEN
1229          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1230          !$ACC PRESENT(tend, tu_m)
1231          DO  i = nxlu, nxr
1232             DO  j = nys, nyn
1233                DO  k = nzb+1, nzt
1234                   tu_m(k,j,i) =    -9.5625_wp * tend(k,j,i)                   &
1235                                   + 5.3125_wp * tu_m(k,j,i)
1236                ENDDO
1237             ENDDO
1238          ENDDO
1239       ENDIF
1240    ENDIF
1241
1242    CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1243
1244!
1245!-- v-velocity component
1246    CALL cpu_log( log_point(6), 'v-equation', 'start' )
1247
1248    !$ACC KERNELS PRESENT(tend)
1249    tend = 0.0_wp
1250    !$ACC END KERNELS
1251    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1252       IF ( ws_scheme_mom )  THEN
1253          CALL advec_v_ws
1254       ELSE
1255          CALL advec_v_pw
1256       END IF
1257    ELSE
1258       CALL advec_v_up
1259    ENDIF
1260    CALL diffusion_v
1261    CALL coriolis( 2 )
1262
1263!
1264!-- Drag by plant canopy
1265    IF ( plant_canopy )  CALL pcm_tendency( 2 )
1266
1267!
1268!-- External pressure gradient
1269    IF ( dp_external )  THEN
1270       DO  i = nxl, nxr
1271          DO  j = nysv, nyn
1272             DO  k = dp_level_ind_b+1, nzt
1273                tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1274             ENDDO
1275          ENDDO
1276       ENDDO
1277    ENDIF
1278
1279!
1280!-- Nudging
1281    IF ( nudging )  CALL nudge( simulated_time, 'v' )
1282
1283!
1284!-- Effect of Stokes drift (in ocean mode only)
1285    IF ( stokes_force )  CALL stokes_drift_terms( 2 )
1286
1287    CALL module_interface_actions( 'v-tendency' )
1288
1289!
1290!-- Prognostic equation for v-velocity component
1291    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1292    !$ACC PRESENT(v, tend, tv_m, v_init, rdf, wall_flags_0) &
1293    !$ACC PRESENT(tsc(2:5)) &
1294    !$ACC PRESENT(v_p)
1295    DO  i = nxl, nxr
1296       DO  j = nysv, nyn
1297          DO  k = nzb+1, nzt
1298             v_p(k,j,i) = v(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1299                                                 tsc(3) * tv_m(k,j,i) )        &
1300                                               - tsc(5) * rdf(k) *             &
1301                                                      ( v(k,j,i) - v_init(k) ) &
1302                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1303                                                BTEST( wall_flags_0(k,j,i), 2 )&
1304                                              )
1305          ENDDO
1306       ENDDO
1307    ENDDO
1308
1309!
1310!-- Add turbulence generated by wave breaking (in ocean mode only)
1311    IF ( wave_breaking  .AND.                                                  &
1312         intermediate_timestep_count == intermediate_timestep_count_max )      &
1313    THEN
1314       CALL wave_breaking_term( 2 )
1315    ENDIF
1316
1317!
1318!-- Calculate tendencies for the next Runge-Kutta step
1319    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1320       IF ( intermediate_timestep_count == 1 )  THEN
1321          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1322          !$ACC PRESENT(tend, tv_m)
1323          DO  i = nxl, nxr
1324             DO  j = nysv, nyn
1325                DO  k = nzb+1, nzt
1326                   tv_m(k,j,i) = tend(k,j,i)
1327                ENDDO
1328             ENDDO
1329          ENDDO
1330       ELSEIF ( intermediate_timestep_count < &
1331                intermediate_timestep_count_max )  THEN
1332          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1333          !$ACC PRESENT(tend, tv_m)
1334          DO  i = nxl, nxr
1335             DO  j = nysv, nyn
1336                DO  k = nzb+1, nzt
1337                   tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1338                                  + 5.3125_wp * tv_m(k,j,i)
1339                ENDDO
1340             ENDDO
1341          ENDDO
1342       ENDIF
1343    ENDIF
1344
1345    CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1346
1347!
1348!-- w-velocity component
1349    CALL cpu_log( log_point(7), 'w-equation', 'start' )
1350
1351    !$ACC KERNELS PRESENT(tend)
1352    tend = 0.0_wp
1353    !$ACC END KERNELS
1354    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1355       IF ( ws_scheme_mom )  THEN
1356          CALL advec_w_ws
1357       ELSE
1358          CALL advec_w_pw
1359       ENDIF
1360    ELSE
1361       CALL advec_w_up
1362    ENDIF
1363    CALL diffusion_w
1364    CALL coriolis( 3 )
1365
1366    IF ( .NOT. neutral )  THEN
1367       IF ( ocean_mode )  THEN
1368          CALL buoyancy( rho_ocean, 3 )
1369       ELSE
1370          IF ( .NOT. humidity )  THEN
1371             CALL buoyancy( pt, 3 )
1372          ELSE
1373             CALL buoyancy( vpt, 3 )
1374          ENDIF
1375       ENDIF
1376    ENDIF
1377
1378!
1379!-- Drag by plant canopy
1380    IF ( plant_canopy )  CALL pcm_tendency( 3 )
1381
1382!
1383!-- Effect of Stokes drift (in ocean mode only)
1384    IF ( stokes_force )  CALL stokes_drift_terms( 3 )
1385
1386    CALL module_interface_actions( 'w-tendency' )
1387
1388!
1389!-- Prognostic equation for w-velocity component
1390    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1391    !$ACC PRESENT(w, tend, tw_m, v_init, rdf, wall_flags_0) &
1392    !$ACC PRESENT(tsc(2:5)) &
1393    !$ACC PRESENT(w_p)
1394    DO  i = nxl, nxr
1395       DO  j = nys, nyn
1396          DO  k = nzb+1, nzt-1
1397             w_p(k,j,i) = w(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1398                                                 tsc(3) * tw_m(k,j,i) )        &
1399                                               - tsc(5) * rdf(k) * w(k,j,i)    &
1400                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1401                                                BTEST( wall_flags_0(k,j,i), 3 )&
1402                                              )
1403          ENDDO
1404       ENDDO
1405    ENDDO
1406
1407!
1408!-- Calculate tendencies for the next Runge-Kutta step
1409    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1410       IF ( intermediate_timestep_count == 1 )  THEN
1411          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1412          !$ACC PRESENT(tend, tw_m)
1413          DO  i = nxl, nxr
1414             DO  j = nys, nyn
1415                DO  k = nzb+1, nzt-1
1416                   tw_m(k,j,i) = tend(k,j,i)
1417                ENDDO
1418             ENDDO
1419          ENDDO
1420       ELSEIF ( intermediate_timestep_count < &
1421                intermediate_timestep_count_max )  THEN
1422          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1423          !$ACC PRESENT(tend, tw_m)
1424          DO  i = nxl, nxr
1425             DO  j = nys, nyn
1426                DO  k = nzb+1, nzt-1
1427                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1428                                  + 5.3125_wp * tw_m(k,j,i)
1429                ENDDO
1430             ENDDO
1431          ENDDO
1432       ENDIF
1433    ENDIF
1434
1435    CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1436
1437
1438!
1439!-- If required, compute prognostic equation for potential temperature
1440    IF ( .NOT. neutral )  THEN
1441
1442       CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1443
1444!
1445!--    pt-tendency terms with communication
1446       sbt = tsc(2)
1447       IF ( scalar_advec == 'bc-scheme' )  THEN
1448
1449          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1450!
1451!--          Bott-Chlond scheme always uses Euler time step. Thus:
1452             sbt = 1.0_wp
1453          ENDIF
1454          tend = 0.0_wp
1455          CALL advec_s_bc( pt, 'pt' )
1456
1457       ENDIF
1458
1459!
1460!--    pt-tendency terms with no communication
1461       IF ( scalar_advec /= 'bc-scheme' )  THEN
1462          !$ACC KERNELS PRESENT(tend)
1463          tend = 0.0_wp
1464          !$ACC END KERNELS
1465          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1466             IF ( ws_scheme_sca )  THEN
1467                CALL advec_s_ws( advc_flags_s, pt, 'pt',                       &
1468                                 bc_dirichlet_l  .OR.  bc_radiation_l,         &
1469                                 bc_dirichlet_n  .OR.  bc_radiation_n,         &
1470                                 bc_dirichlet_r  .OR.  bc_radiation_r,         &
1471                                 bc_dirichlet_s  .OR.  bc_radiation_s )
1472             ELSE
1473                CALL advec_s_pw( pt )
1474             ENDIF
1475          ELSE
1476             CALL advec_s_up( pt )
1477          ENDIF
1478       ENDIF
1479
1480       CALL diffusion_s( pt,                                                   &
1481                         surf_def_h(0)%shf, surf_def_h(1)%shf,                 &
1482                         surf_def_h(2)%shf,                                    &
1483                         surf_lsm_h%shf,    surf_usm_h%shf,                    &
1484                         surf_def_v(0)%shf, surf_def_v(1)%shf,                 &
1485                         surf_def_v(2)%shf, surf_def_v(3)%shf,                 &
1486                         surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,                 &
1487                         surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,                 &
1488                         surf_usm_v(0)%shf, surf_usm_v(1)%shf,                 &
1489                         surf_usm_v(2)%shf, surf_usm_v(3)%shf )
1490
1491!
1492!--    Consideration of heat sources within the plant canopy
1493       IF ( plant_canopy  .AND.                                          &
1494            (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
1495          CALL pcm_tendency( 4 )
1496       ENDIF
1497
1498!
1499!--    Large scale advection
1500       IF ( large_scale_forcing )  THEN
1501          CALL ls_advec( simulated_time, 'pt' )
1502       ENDIF
1503
1504!
1505!--    Nudging
1506       IF ( nudging )  CALL nudge( simulated_time, 'pt' )
1507
1508!
1509!--    If required compute influence of large-scale subsidence/ascent
1510       IF ( large_scale_subsidence  .AND.                                      &
1511            .NOT. use_subsidence_tendencies )  THEN
1512          CALL subsidence( tend, pt, pt_init, 2 )
1513       ENDIF
1514
1515#if defined( __rrtmg )
1516!
1517!--    If required, add tendency due to radiative heating/cooling
1518       IF ( radiation  .AND.                                                   &
1519            simulated_time > skip_time_do_radiation )  THEN
1520            CALL radiation_tendency ( tend )
1521       ENDIF
1522#endif
1523
1524       CALL module_interface_actions( 'pt-tendency' )
1525
1526!
1527!--    Prognostic equation for potential temperature
1528       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1529       !$ACC PRESENT(pt, tend, tpt_m, wall_flags_0) &
1530       !$ACC PRESENT(pt_init, rdf_sc, ptdf_x, ptdf_y) &
1531       !$ACC PRESENT(tsc(3:5)) &
1532       !$ACC PRESENT(pt_p)
1533       DO  i = nxl, nxr
1534          DO  j = nys, nyn
1535             DO  k = nzb+1, nzt
1536                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +      &
1537                                                   tsc(3) * tpt_m(k,j,i) )     &
1538                                                 - tsc(5) *                    &
1539                                                   ( pt(k,j,i) - pt_init(k) ) *&
1540                                          ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )&
1541                                          )                                    &
1542                                   * MERGE( 1.0_wp, 0.0_wp,                    &
1543                                             BTEST( wall_flags_0(k,j,i), 0 )   &
1544                                          )
1545             ENDDO
1546          ENDDO
1547       ENDDO
1548!
1549!--    Calculate tendencies for the next Runge-Kutta step
1550       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1551          IF ( intermediate_timestep_count == 1 )  THEN
1552             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1553             !$ACC PRESENT(tend, tpt_m)
1554             DO  i = nxl, nxr
1555                DO  j = nys, nyn
1556                   DO  k = nzb+1, nzt
1557                      tpt_m(k,j,i) = tend(k,j,i)
1558                   ENDDO
1559                ENDDO
1560             ENDDO
1561          ELSEIF ( intermediate_timestep_count < &
1562                   intermediate_timestep_count_max )  THEN
1563             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1564             !$ACC PRESENT(tend, tpt_m)
1565             DO  i = nxl, nxr
1566                DO  j = nys, nyn
1567                   DO  k = nzb+1, nzt
1568                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
1569                                        5.3125_wp * tpt_m(k,j,i)
1570                   ENDDO
1571                ENDDO
1572             ENDDO
1573          ENDIF
1574       ENDIF
1575
1576       CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1577
1578    ENDIF
1579
1580!
1581!-- If required, compute prognostic equation for total water content
1582    IF ( humidity )  THEN
1583
1584       CALL cpu_log( log_point(29), 'q-equation', 'start' )
1585
1586!
1587!--    Scalar/q-tendency terms with communication
1588       sbt = tsc(2)
1589       IF ( scalar_advec == 'bc-scheme' )  THEN
1590
1591          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1592!
1593!--          Bott-Chlond scheme always uses Euler time step. Thus:
1594             sbt = 1.0_wp
1595          ENDIF
1596          tend = 0.0_wp
1597          CALL advec_s_bc( q, 'q' )
1598
1599       ENDIF
1600
1601!
1602!--    Scalar/q-tendency terms with no communication
1603       IF ( scalar_advec /= 'bc-scheme' )  THEN
1604          tend = 0.0_wp
1605          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1606             IF ( ws_scheme_sca )  THEN
1607                CALL advec_s_ws( advc_flags_s, q, 'q',                         &
1608                                 bc_dirichlet_l  .OR.  bc_radiation_l,         &
1609                                 bc_dirichlet_n  .OR.  bc_radiation_n,         &
1610                                 bc_dirichlet_r  .OR.  bc_radiation_r,         &
1611                                 bc_dirichlet_s  .OR.  bc_radiation_s )
1612             ELSE
1613                CALL advec_s_pw( q )
1614             ENDIF
1615          ELSE
1616             CALL advec_s_up( q )
1617          ENDIF
1618       ENDIF
1619
1620       CALL diffusion_s( q,                                                    &
1621                         surf_def_h(0)%qsws, surf_def_h(1)%qsws,               &
1622                         surf_def_h(2)%qsws,                                   &
1623                         surf_lsm_h%qsws,    surf_usm_h%qsws,                  &
1624                         surf_def_v(0)%qsws, surf_def_v(1)%qsws,               &
1625                         surf_def_v(2)%qsws, surf_def_v(3)%qsws,               &
1626                         surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,               &
1627                         surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,               &
1628                         surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,               &
1629                         surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
1630
1631!
1632!--    Sink or source of humidity due to canopy elements
1633       IF ( plant_canopy ) CALL pcm_tendency( 5 )
1634
1635!
1636!--    Large scale advection
1637       IF ( large_scale_forcing )  THEN
1638          CALL ls_advec( simulated_time, 'q' )
1639       ENDIF
1640
1641!
1642!--    Nudging
1643       IF ( nudging )  CALL nudge( simulated_time, 'q' )
1644
1645!
1646!--    If required compute influence of large-scale subsidence/ascent
1647       IF ( large_scale_subsidence  .AND.                                      &
1648            .NOT. use_subsidence_tendencies )  THEN
1649         CALL subsidence( tend, q, q_init, 3 )
1650       ENDIF
1651
1652       CALL module_interface_actions( 'q-tendency' )
1653
1654!
1655!--    Prognostic equation for total water content
1656       DO  i = nxl, nxr
1657          DO  j = nys, nyn
1658             DO  k = nzb+1, nzt
1659                q_p(k,j,i) = q(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
1660                                                 tsc(3) * tq_m(k,j,i) )        &
1661                                               - tsc(5) * rdf_sc(k) *          &
1662                                                      ( q(k,j,i) - q_init(k) ) &
1663                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
1664                                                   BTEST( wall_flags_0(k,j,i), 0 ) &
1665                                                 )
1666                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
1667             ENDDO
1668          ENDDO
1669       ENDDO
1670
1671!
1672!--    Calculate tendencies for the next Runge-Kutta step
1673       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1674          IF ( intermediate_timestep_count == 1 )  THEN
1675             DO  i = nxl, nxr
1676                DO  j = nys, nyn
1677                   DO  k = nzb+1, nzt
1678                      tq_m(k,j,i) = tend(k,j,i)
1679                   ENDDO
1680                ENDDO
1681             ENDDO
1682          ELSEIF ( intermediate_timestep_count < &
1683                   intermediate_timestep_count_max )  THEN
1684             DO  i = nxl, nxr
1685                DO  j = nys, nyn
1686                   DO  k = nzb+1, nzt
1687                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
1688                                     + 5.3125_wp * tq_m(k,j,i)
1689                   ENDDO
1690                ENDDO
1691             ENDDO
1692          ENDIF
1693       ENDIF
1694
1695       CALL cpu_log( log_point(29), 'q-equation', 'stop' )
1696
1697    ENDIF
1698!
1699!-- If required, compute prognostic equation for scalar
1700    IF ( passive_scalar )  THEN
1701
1702       CALL cpu_log( log_point(66), 's-equation', 'start' )
1703
1704!
1705!--    Scalar/q-tendency terms with communication
1706       sbt = tsc(2)
1707       IF ( scalar_advec == 'bc-scheme' )  THEN
1708
1709          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1710!
1711!--          Bott-Chlond scheme always uses Euler time step. Thus:
1712             sbt = 1.0_wp
1713          ENDIF
1714          tend = 0.0_wp
1715          CALL advec_s_bc( s, 's' )
1716
1717       ENDIF
1718
1719!
1720!--    Scalar/q-tendency terms with no communication
1721       IF ( scalar_advec /= 'bc-scheme' )  THEN
1722          tend = 0.0_wp
1723          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1724             IF ( ws_scheme_sca )  THEN
1725                CALL advec_s_ws( advc_flags_s, s, 's',                         &
1726                                 bc_dirichlet_l  .OR.  bc_radiation_l,         &
1727                                 bc_dirichlet_n  .OR.  bc_radiation_n,         &
1728                                 bc_dirichlet_r  .OR.  bc_radiation_r,         &
1729                                 bc_dirichlet_s  .OR.  bc_radiation_s )
1730             ELSE
1731                CALL advec_s_pw( s )
1732             ENDIF
1733          ELSE
1734             CALL advec_s_up( s )
1735          ENDIF
1736       ENDIF
1737
1738       CALL diffusion_s( s,                                                    &
1739                         surf_def_h(0)%ssws, surf_def_h(1)%ssws,               &
1740                         surf_def_h(2)%ssws,                                   &
1741                         surf_lsm_h%ssws,    surf_usm_h%ssws,                  &
1742                         surf_def_v(0)%ssws, surf_def_v(1)%ssws,               &
1743                         surf_def_v(2)%ssws, surf_def_v(3)%ssws,               &
1744                         surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,               &
1745                         surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,               &
1746                         surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,               &
1747                         surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1748
1749!
1750!--    Sink or source of humidity due to canopy elements
1751       IF ( plant_canopy ) CALL pcm_tendency( 7 )
1752
1753!
1754!--    Large scale advection. Not implemented for scalars so far.
1755!        IF ( large_scale_forcing )  THEN
1756!           CALL ls_advec( simulated_time, 'q' )
1757!        ENDIF
1758
1759!
1760!--    Nudging. Not implemented for scalars so far.
1761!        IF ( nudging )  CALL nudge( simulated_time, 'q' )
1762
1763!
1764!--    If required compute influence of large-scale subsidence/ascent.
1765!--    Not implemented for scalars so far.
1766       IF ( large_scale_subsidence  .AND.                                      &
1767            .NOT. use_subsidence_tendencies  .AND.                             &
1768            .NOT. large_scale_forcing )  THEN
1769         CALL subsidence( tend, s, s_init, 3 )
1770       ENDIF
1771
1772       CALL module_interface_actions( 's-tendency' )
1773
1774!
1775!--    Prognostic equation for total water content
1776       DO  i = nxl, nxr
1777          DO  j = nys, nyn
1778             DO  k = nzb+1, nzt
1779                s_p(k,j,i) = s(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
1780                                                 tsc(3) * ts_m(k,j,i) )        &
1781                                               - tsc(5) * rdf_sc(k) *          &
1782                                                    ( s(k,j,i) - s_init(k) )   &
1783                                        )                                      &
1784                                   * MERGE( 1.0_wp, 0.0_wp,                    &
1785                                             BTEST( wall_flags_0(k,j,i), 0 )   &
1786                                          )
1787                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
1788             ENDDO
1789          ENDDO
1790       ENDDO
1791
1792!
1793!--    Calculate tendencies for the next Runge-Kutta step
1794       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1795          IF ( intermediate_timestep_count == 1 )  THEN
1796             DO  i = nxl, nxr
1797                DO  j = nys, nyn
1798                   DO  k = nzb+1, nzt
1799                      ts_m(k,j,i) = tend(k,j,i)
1800                   ENDDO
1801                ENDDO
1802             ENDDO
1803          ELSEIF ( intermediate_timestep_count < &
1804                   intermediate_timestep_count_max )  THEN
1805             DO  i = nxl, nxr
1806                DO  j = nys, nyn
1807                   DO  k = nzb+1, nzt
1808                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
1809                                     + 5.3125_wp * ts_m(k,j,i)
1810                   ENDDO
1811                ENDDO
1812             ENDDO
1813          ENDIF
1814       ENDIF
1815
1816       CALL cpu_log( log_point(66), 's-equation', 'stop' )
1817
1818    ENDIF
1819!
1820!-- Calculate prognostic equations for all other modules
1821    CALL module_interface_prognostic_equations()
1822
1823    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_vector', 'end' )
1824
1825 END SUBROUTINE prognostic_equations_vector
1826
1827
1828 END MODULE prognostic_equations_mod
Note: See TracBrowser for help on using the repository browser.