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

Last change on this file since 4723 was 4717, checked in by pavelkrc, 4 years ago

Fixes and optimizations of OpenMP parallelization, formatting of OpenMP directives

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