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

Last change on this file since 4360 was 4360, checked in by suehring, 16 months ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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