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

Last change on this file since 4079 was 4079, checked in by suehring, 5 years ago

Implementation of a monotonic flux limiter for vertical advection term in Wicker-Skamarock scheme. The flux limiter is currently only applied for passive scalars (passive scalar, chemical species, aerosols) within the region up to the highest topography, in order to avoid the built-up of large concentrations within poorly resolved cavities in urban environments. To enable the limiter monotonic_limiter_z = .T. must be set. Note, the limiter is currently only implemented for the cache-optimized version of advec_ws. Further changes in offline nesting: Set boundary condition for w at nzt+1 at all lateral boundaries (even though these won't enter the numerical solution), in order to avoid high vertical velocities in the run-control file which might built-up due to the mass-conservation; bugfix in offline nesting for chemical species

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