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

Last change on this file since 4719 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
Line 
1!> @file prognostic_equations.f90
2!--------------------------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the 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.
8!
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.
12!
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/>.
15!
16! Copyright 1997-2020 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: prognostic_equations.f90 4717 2020-09-30 22:27:40Z pavelkrc $
27! Fixes and optimizations of OpenMP parallelization, formatting of OpenMP
28! directives (J. Resler)
29!
30! 4671 2020-09-09 20:27:58Z pavelkrc
31! Implementation of downward facing USM and LSM surfaces
32!
33! 4649 2020-08-25 12:11:17Z raasch
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!
39! 4360 2020-01-07 11:25:50Z suehring
40! Introduction of wall_flags_total_0, which currently sets bits based on static topography
41! information used in wall_flags_static_0
42!
43! 4329 2019-12-10 15:46:36Z motisi
44! Renamed wall_flags_0 to wall_flags_static_0
45!
46! 4182 2019-08-22 15:20:23Z scharf
47! Corrected "Former revisions" section
48!
49! 4110 2019-07-22 17:05:21Z suehring
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!
54! 4109 2019-07-22 17:00:34Z suehring
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!
59! 4048 2019-06-21 21:00:21Z knoop
60! Moved tcm_prognostic_equations to module_interface
61!
62! 3987 2019-05-22 09:52:13Z kanani
63! Introduce alternative switch for debug output during timestepping
64!
65! 3956 2019-05-07 12:32:52Z monakurppa
66! Removed salsa calls.
67!
68! 3931 2019-04-24 16:34:28Z schwenkel
69! Correct/complete module_interface introduction for chemistry model
70!
71! 3899 2019-04-16 14:05:27Z monakurppa
72! Corrections in the OpenMP version of salsa
73!
74! 3887 2019 -04-12 08:47:41Z schwenkel
75! Implicit Bugfix for chemistry model, loop for non_transport_physics over ghost points is avoided.
76! Instead introducing module_interface_exchange_horiz.
77!
78! 3885 2019-04-11 11:29:34Z kanani
79! Changes related to global restructuring of location messages and introduction of additional debug
80! messages
81!
82! 3881 2019-04-10 09:31:22Z suehring
83! Bugfix in OpenMP directive
84!
85! 3880 2019-04-08 21:43:02Z knoop
86! Moved wtm_tendencies to module_interface_actions
87!
88! 3874 2019-04-08 16:53:48Z knoop
89! Added non_transport_physics module interfaces and moved bcm code into it
90!
91! 3872 2019-04-08 15:03:06Z knoop
92! Moving prognostic equations of bcm into bulk_cloud_model_mod
93!
94! 3864 2019-04-05 09:01:56Z monakurppa
95! Modifications made for salsa:
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!
100! 3840 2019-03-29 10:35:52Z knoop
101! Added USE chem_gasphase_mod for nspec, nspec and spc_names
102!
103! 3820 2019-03-27 11:53:41Z forkel
104! Renamed do_depo to deposition_dry (ecc)
105!
106! 3797 2019-03-15 11:15:38Z forkel
107! Call chem_integegrate in OpenMP loop (ketelsen)
108!
109!
110! 3771 2019-02-28 12:19:33Z raasch
111! Preprocessor directivs fro rrtmg added
112!
113! 3761 2019-02-25 15:31:42Z raasch
114! Unused variable removed
115!
116! 3719 2019-02-06 13:10:18Z kanani
117! Cleaned up chemistry cpu measurements
118!
119! 3684 2019-01-20 20:20:58Z knoop
120! OpenACC port for SPEC
121!
122! Revision 1.1  2000/04/13 14:56:27  schroeter
123! Initial revision
124!
125!--------------------------------------------------------------------------------------------------!
126! Description:
127! ------------
128!> Solving the prognostic equations.
129!--------------------------------------------------------------------------------------------------!
130 MODULE prognostic_equations_mod
131
132    USE advec_s_bc_mod,                                                                            &
133        ONLY:  advec_s_bc
134
135    USE advec_s_pw_mod,                                                                            &
136        ONLY:  advec_s_pw
137
138    USE advec_s_up_mod,                                                                            &
139        ONLY:  advec_s_up
140
141    USE advec_u_pw_mod,                                                                            &
142        ONLY:  advec_u_pw
143
144    USE advec_u_up_mod,                                                                            &
145        ONLY:  advec_u_up
146
147    USE advec_v_pw_mod,                                                                            &
148        ONLY:  advec_v_pw
149
150    USE advec_v_up_mod,                                                                            &
151        ONLY:  advec_v_up
152
153    USE advec_w_pw_mod,                                                                            &
154        ONLY:  advec_w_pw
155
156    USE advec_w_up_mod,                                                                            &
157        ONLY:  advec_w_up
158
159    USE advec_ws,                                                                                  &
160        ONLY:  advec_s_ws,                                                                         &
161               advec_u_ws,                                                                         &
162               advec_v_ws,                                                                         &
163               advec_w_ws
164
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
223
224    USE buoyancy_mod,                                                                              &
225        ONLY:  buoyancy
226
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
270
271
272
273
274
275    USE coriolis_mod,                                                                              &
276        ONLY:  coriolis
277
278    USE cpulog,                                                                                    &
279        ONLY:  cpu_log,                                                                            &
280               log_point,                                                                          &
281               log_point_s
282
283    USE diffusion_s_mod,                                                                           &
284        ONLY:  diffusion_s
285
286    USE diffusion_u_mod,                                                                           &
287        ONLY:  diffusion_u
288
289    USE diffusion_v_mod,                                                                           &
290        ONLY:  diffusion_v
291
292    USE diffusion_w_mod,                                                                           &
293        ONLY:  diffusion_w
294
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
311
312    USE kinds
313
314    USE lsf_nudging_mod,                                                                           &
315        ONLY:  ls_advec,                                                                           &
316               nudge
317
318    USE module_interface,                                                                          &
319        ONLY:  module_interface_actions,                                                           &
320               module_interface_exchange_horiz,                                                    &
321               module_interface_non_advective_processes,                                           &
322               module_interface_prognostic_equations
323
324    USE ocean_mod,                                                                                 &
325        ONLY:  stokes_drift_terms,                                                                 &
326               stokes_force,                                                                       &
327               wave_breaking,                                                                      &
328               wave_breaking_term
329
330    USE plant_canopy_model_mod,                                                                    &
331        ONLY:  cthf,                                                                               &
332               pcm_tendency
333
334#if defined( __rrtmg )
335    USE radiation_model_mod,                                                                       &
336        ONLY:  radiation,                                                                          &
337               radiation_tendency,                                                                 &
338               skip_time_do_radiation
339#endif
340
341    USE statistics,                                                                                &
342        ONLY:  hom
343
344    USE subsidence_mod,                                                                            &
345        ONLY:  subsidence
346
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,                                                                        &
353                surf_usm_v
354
355    IMPLICIT NONE
356
357    PRIVATE
358    PUBLIC prognostic_equations_cache, prognostic_equations_vector
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
372!--------------------------------------------------------------------------------------------------!
373! Description:
374! ------------
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.
377!>
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.
380!>
381!> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
382!--------------------------------------------------------------------------------------------------!
383
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                   !<
391!$  INTEGER(iwp) ::  omp_get_thread_num  !<
392    INTEGER(iwp) ::  tn = 0              !<
393
394    LOGICAL      ::  loop_start          !<
395
396
397    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_cache', 'start' )
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
402    !$OMP PARALLEL DO PRIVATE (i,j) SCHEDULE( STATIC )
403    DO  i = nxl, nxr
404       DO  j = nys, nyn
405!
406!--       Calculate non advective processes for all other modules
407          CALL module_interface_non_advective_processes( i, j )
408       ENDDO
409    ENDDO
410!
411!-- Module Inferface for exchange horiz after non_advective_processes but before advection.
412!-- Therefore, non_advective_processes must not run for ghost points.
413    CALL module_interface_exchange_horiz()
414!
415!-- Loop over all prognostic equations
416    !$OMP PARALLEL PRIVATE (i,i_omp_start,j,k,loop_start,tn)
417
418    !$  tn = omp_get_thread_num()
419    loop_start = .TRUE.
420
421    !$OMP DO
422    DO  i = nxl, nxr
423
424!
425!--    Store the first loop index. It differs for each thread and is required later in advec_ws
426       IF ( loop_start )  THEN
427          loop_start  = .FALSE.
428          i_omp_start = i
429       ENDIF
430
431       DO  j = nys, nyn
432!
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.
436          IF ( i >= nxlu )  THEN
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 )
442                ELSE
443                   CALL advec_u_pw( i, j )
444                ENDIF
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
470!
471!--          Effect of Stokes drift (in ocean mode only)
472             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 1 )
473
474             CALL module_interface_actions( i, j, 'u-tendency' )
475!
476!--          Prognostic equation for u-velocity component
477             DO  k = nzb+1, nzt
478
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 ) )
482             ENDDO
483
484!
485!--          Add turbulence generated by wave breaking (in ocean mode only)
486             IF ( wave_breaking  .AND.                                                             &
487                  intermediate_timestep_count == intermediate_timestep_count_max )  THEN
488                CALL wave_breaking_term( i, j, 1 )
489             ENDIF
490
491!
492!--          Calculate tendencies for the next Runge-Kutta step
493             IF ( timestep_scheme(1:5) == 'runge' )  THEN
494                IF ( intermediate_timestep_count == 1 )  THEN
495                   DO  k = nzb+1, nzt
496                      tu_m(k,j,i) = tend(k,j,i)
497                   ENDDO
498                ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
499                   DO  k = nzb+1, nzt
500                      tu_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tu_m(k,j,i)
501                   ENDDO
502                ENDIF
503             ENDIF
504
505          ENDIF
506!
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.
510          IF ( j >= nysv )  THEN
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 )
516                ELSE
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
527             IF ( plant_canopy )  CALL pcm_tendency( i, j, 2 )
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
541!
542!--          Effect of Stokes drift (in ocean mode only)
543             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 2 )
544
545             CALL module_interface_actions( i, j, 'v-tendency' )
546!
547!--          Prognostic equation for v-velocity component
548             DO  k = nzb+1, nzt
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 ) )
552             ENDDO
553
554!
555!--          Add turbulence generated by wave breaking (in ocean mode only)
556             IF ( wave_breaking  .AND.                                                             &
557                  intermediate_timestep_count == intermediate_timestep_count_max )  THEN
558                CALL wave_breaking_term( i, j, 2 )
559             ENDIF
560
561!
562!--          Calculate tendencies for the next Runge-Kutta step
563             IF ( timestep_scheme(1:5) == 'runge' )  THEN
564                IF ( intermediate_timestep_count == 1 )  THEN
565                   DO  k = nzb+1, nzt
566                      tv_m(k,j,i) = tend(k,j,i)
567                   ENDDO
568                ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
569                   DO  k = nzb+1, nzt
570                      tv_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tv_m(k,j,i)
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 )
583             ELSE
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
593             IF ( ocean_mode )  THEN
594                CALL buoyancy( i, j, rho_ocean, 3 )
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
608!
609!--       Effect of Stokes drift (in ocean mode only)
610          IF ( stokes_force )  CALL stokes_drift_terms( i, j, 3 )
611
612          CALL module_interface_actions( i, j, 'w-tendency' )
613!
614!--       Prognostic equation for w-velocity component
615          DO  k = nzb+1, nzt-1
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 ) )
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
625                DO  k = nzb+1, nzt-1
626                   tw_m(k,j,i) = tend(k,j,i)
627                ENDDO
628             ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
629                DO  k = nzb+1, nzt-1
630                   tw_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tw_m(k,j,i)
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
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,                       &
648                                       bc_dirichlet_s  .OR.  bc_radiation_s )
649                   ELSE
650                      CALL advec_s_pw( i, j, pt )
651                   ENDIF
652             ELSE
653                CALL advec_s_up( i, j, pt )
654             ENDIF
655             CALL diffusion_s( i, j, pt, surf_def_h(0)%shf, surf_def_h(1)%shf, surf_def_h(2)%shf,  &
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 )
662
663!
664!--          Consideration of heat sources within the plant canopy
665             IF ( plant_canopy  .AND.  (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  &
666             THEN
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' )
674             ENDIF
675
676!
677!--          Nudging
678             IF ( nudging )  CALL nudge( i, j, simulated_time, 'pt' )
679
680!
681!--          If required, compute effect of large-scale subsidence/ascent
682             IF ( large_scale_subsidence  .AND.  .NOT. use_subsidence_tendencies )  THEN
683                CALL subsidence( i, j, tend, pt, pt_init, 2 )
684             ENDIF
685
686#if defined( __rrtmg )
687!
688!--          If required, add tendency due to radiative heating/cooling
689             IF ( radiation  .AND.  simulated_time > skip_time_do_radiation )  THEN
690                CALL radiation_tendency ( i, j, tend )
691             ENDIF
692#endif
693
694             CALL module_interface_actions( i, j, 'pt-tendency' )
695!
696!--          Prognostic equation for potential temperature
697             DO  k = nzb+1, nzt
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 ) )
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
709                   DO  k = nzb+1, nzt
710                      tpt_m(k,j,i) = tend(k,j,i)
711                   ENDDO
712                ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
713                   DO  k = nzb+1, nzt
714                      tpt_m(k,j,i) = - 9.5625_wp * tend(k,j,i) + 5.3125_wp * tpt_m(k,j,i)
715                   ENDDO
716                ENDIF
717             ENDIF
718
719          ENDIF
720
721!
722!--       If required, compute prognostic equation for total water content
723          IF ( humidity )  THEN
724
725!
726!--          Tendency-terms for total water content / scalar
727             tend(:,j,i) = 0.0_wp
728             IF ( timestep_scheme(1:5) == 'runge' )  THEN
729                IF ( ws_scheme_sca )  THEN
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,                          &
735                                    bc_dirichlet_s  .OR.  bc_radiation_s )
736                ELSE
737                   CALL advec_s_pw( i, j, q )
738                ENDIF
739             ELSE
740                CALL advec_s_up( i, j, q )
741             ENDIF
742             CALL diffusion_s( i, j, q, surf_def_h(0)%qsws, surf_def_h(1)%qsws,                    &
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,                             &
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 )
749
750!
751!--          Sink or source of humidity due to canopy elements
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
762             IF ( nudging )  CALL nudge( i, j, simulated_time, 'q' )
763
764!
765!--          If required compute influence of large-scale subsidence/ascent
766             IF ( large_scale_subsidence  .AND.  .NOT. use_subsidence_tendencies )  THEN
767                CALL subsidence( i, j, tend, q, q_init, 3 )
768             ENDIF
769
770             CALL module_interface_actions( i, j, 'q-tendency' )
771
772!
773!--          Prognostic equation for total water content / scalar
774             DO  k = nzb+1, nzt
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 ) )
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
786                   DO  k = nzb+1, nzt
787                      tq_m(k,j,i) = tend(k,j,i)
788                   ENDDO
789                ELSEIF ( intermediate_timestep_count < &
790                         intermediate_timestep_count_max )  THEN
791                   DO  k = nzb+1, nzt
792                      tq_m(k,j,i) = - 9.5625_wp * tend(k,j,i) + 5.3125_wp * tq_m(k,j,i)
793                   ENDDO
794                ENDIF
795             ENDIF
796
797          ENDIF
798
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
805             IF ( timestep_scheme(1:5) == 'runge' )  THEN
806                IF ( ws_scheme_sca )  THEN
807!
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,                          &
815                                    monotonic_limiter_z )
816                ELSE
817                   CALL advec_s_pw( i, j, s )
818                ENDIF
819             ELSE
820                CALL advec_s_up( i, j, s )
821             ENDIF
822             CALL diffusion_s( i, j, s, surf_def_h(0)%ssws, surf_def_h(1)%ssws,                    &
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,                             &
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 )
829
830!
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
842!              IF ( nudging )  CALL nudge( i, j, simulated_time, 's' )
843
844!
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.            &
849                  .NOT. large_scale_forcing )  THEN
850                CALL subsidence( i, j, tend, s, s_init, 3 )
851             ENDIF
852
853             CALL module_interface_actions( i, j, 's-tendency' )
854
855!
856!--          Prognostic equation for scalar
857             DO  k = nzb+1, nzt
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 ) )
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
869                   DO  k = nzb+1, nzt
870                      ts_m(k,j,i) = tend(k,j,i)
871                   ENDDO
872                ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
873                   DO  k = nzb+1, nzt
874                      ts_m(k,j,i) = - 9.5625_wp * tend(k,j,i) + 5.3125_wp * ts_m(k,j,i)
875                   ENDDO
876                ENDIF
877             ENDIF
878
879          ENDIF
880!
881!--       Calculate prognostic equations for all other modules
882          CALL module_interface_prognostic_equations( i, j, i_omp_start, tn )
883
884       ENDDO  ! loop over j
885    ENDDO  ! loop over i
886    !$OMP END PARALLEL
887
888
889    CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
890
891    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_cache', 'end' )
892
893 END SUBROUTINE prognostic_equations_cache
894
895
896!--------------------------------------------------------------------------------------------------!
897! Description:
898! ------------
899!> Version for vector machines
900!--------------------------------------------------------------------------------------------------!
901
902 SUBROUTINE prognostic_equations_vector
903
904
905    INTEGER(iwp) ::  i  !<
906    INTEGER(iwp) ::  j  !<
907    INTEGER(iwp) ::  k  !<
908
909    REAL(wp) ::  sbt  !<
910
911
912    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_vector', 'start' )
913!
914!-- Calculate non advective processes for all other modules
915    CALL module_interface_non_advective_processes
916!
917!-- Module Inferface for exchange horiz after non_advective_processes but before advection.
918!-- Therefore, non_advective_processes must not run for ghost points.
919    CALL module_interface_exchange_horiz()
920!
921!-- u-velocity component
922    CALL cpu_log( log_point(5), 'u-equation', 'start' )
923
924    !$ACC KERNELS PRESENT(tend)
925    tend = 0.0_wp
926    !$ACC END KERNELS
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
962!
963!-- Effect of Stokes drift (in ocean mode only)
964    IF ( stokes_force )  CALL stokes_drift_terms( 1 )
965
966    CALL module_interface_actions( 'u-tendency' )
967
968!
969!-- Prognostic equation for u-velocity component
970    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
971    !$ACC PRESENT(u, tend, tu_m, u_init, rdf, wall_flags_total_0) &
972    !$ACC PRESENT(tsc(2:5)) &
973    !$ACC PRESENT(u_p)
974    DO  i = nxlu, nxr
975       DO  j = nys, nyn
976!
977!--       Following directive is required to vectorize on Intel19
978          !DIR$ IVDEP
979          DO  k = nzb+1, nzt
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 ) )
984          ENDDO
985       ENDDO
986    ENDDO
987
988!
989!-- Add turbulence generated by wave breaking (in ocean mode only)
990    IF ( wave_breaking .AND. intermediate_timestep_count == intermediate_timestep_count_max )  THEN
991       CALL wave_breaking_term( 1 )
992    ENDIF
993
994!
995!-- Calculate tendencies for the next Runge-Kutta step
996    IF ( timestep_scheme(1:5) == 'runge' )  THEN
997       IF ( intermediate_timestep_count == 1 )  THEN
998          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
999          !$ACC PRESENT(tend, tu_m)
1000          DO  i = nxlu, nxr
1001             DO  j = nys, nyn
1002                DO  k = nzb+1, nzt
1003                   tu_m(k,j,i) = tend(k,j,i)
1004                ENDDO
1005             ENDDO
1006          ENDDO
1007       ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
1008          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1009          !$ACC PRESENT(tend, tu_m)
1010          DO  i = nxlu, nxr
1011             DO  j = nys, nyn
1012                DO  k = nzb+1, nzt
1013                   tu_m(k,j,i) = - 9.5625_wp * tend(k,j,i) + 5.3125_wp * tu_m(k,j,i)
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
1026    !$ACC KERNELS PRESENT(tend)
1027    tend = 0.0_wp
1028    !$ACC END KERNELS
1029    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1030       IF ( ws_scheme_mom )  THEN
1031          CALL advec_v_ws
1032       ELSE
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
1061!
1062!-- Effect of Stokes drift (in ocean mode only)
1063    IF ( stokes_force )  CALL stokes_drift_terms( 2 )
1064
1065    CALL module_interface_actions( 'v-tendency' )
1066
1067!
1068!-- Prognostic equation for v-velocity component
1069    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1070    !$ACC PRESENT(v, tend, tv_m, v_init, rdf, wall_flags_total_0) &
1071    !$ACC PRESENT(tsc(2:5)) &
1072    !$ACC PRESENT(v_p)
1073    DO  i = nxl, nxr
1074       DO  j = nysv, nyn
1075!
1076!--       Following directive is required to vectorize on Intel19
1077          !DIR$ IVDEP
1078          DO  k = nzb+1, nzt
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 ) )
1083          ENDDO
1084       ENDDO
1085    ENDDO
1086
1087!
1088!-- Add turbulence generated by wave breaking (in ocean mode only)
1089    IF ( wave_breaking .AND. intermediate_timestep_count == intermediate_timestep_count_max )  THEN
1090       CALL wave_breaking_term( 2 )
1091    ENDIF
1092
1093!
1094!-- Calculate tendencies for the next Runge-Kutta step
1095    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1096       IF ( intermediate_timestep_count == 1 )  THEN
1097          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1098          !$ACC PRESENT(tend, tv_m)
1099          DO  i = nxl, nxr
1100             DO  j = nysv, nyn
1101                DO  k = nzb+1, nzt
1102                   tv_m(k,j,i) = tend(k,j,i)
1103                ENDDO
1104             ENDDO
1105          ENDDO
1106       ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
1107          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1108          !$ACC PRESENT(tend, tv_m)
1109          DO  i = nxl, nxr
1110             DO  j = nysv, nyn
1111                DO  k = nzb+1, nzt
1112                   tv_m(k,j,i) = - 9.5625_wp * tend(k,j,i) + 5.3125_wp * tv_m(k,j,i)
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
1125    !$ACC KERNELS PRESENT(tend)
1126    tend = 0.0_wp
1127    !$ACC END KERNELS
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
1141       IF ( ocean_mode )  THEN
1142          CALL buoyancy( rho_ocean, 3 )
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
1156!
1157!-- Effect of Stokes drift (in ocean mode only)
1158    IF ( stokes_force )  CALL stokes_drift_terms( 3 )
1159
1160    CALL module_interface_actions( 'w-tendency' )
1161
1162!
1163!-- Prognostic equation for w-velocity component
1164    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1165    !$ACC PRESENT(w, tend, tw_m, v_init, rdf, wall_flags_total_0) &
1166    !$ACC PRESENT(tsc(2:5)) &
1167    !$ACC PRESENT(w_p)
1168    DO  i = nxl, nxr
1169       DO  j = nys, nyn
1170!
1171!--       Following directive is required to vectorize on Intel19
1172          !DIR$ IVDEP
1173          DO  k = nzb+1, nzt-1
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 ) )
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
1186          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1187          !$ACC PRESENT(tend, tw_m)
1188          DO  i = nxl, nxr
1189             DO  j = nys, nyn
1190                DO  k = nzb+1, nzt-1
1191                   tw_m(k,j,i) = tend(k,j,i)
1192                ENDDO
1193             ENDDO
1194          ENDDO
1195       ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
1196          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1197          !$ACC PRESENT(tend, tw_m)
1198          DO  i = nxl, nxr
1199             DO  j = nys, nyn
1200                DO  k = nzb+1, nzt-1
1201                   tw_m(k,j,i) = - 9.5625_wp * tend(k,j,i) + 5.3125_wp * tw_m(k,j,i)
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
1235          !$ACC KERNELS PRESENT(tend)
1236          tend = 0.0_wp
1237          !$ACC END KERNELS
1238          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1239             IF ( ws_scheme_sca )  THEN
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,                             &
1244                                 bc_dirichlet_s  .OR.  bc_radiation_s )
1245             ELSE
1246                CALL advec_s_pw( pt )
1247             ENDIF
1248          ELSE
1249             CALL advec_s_up( pt )
1250          ENDIF
1251       ENDIF
1252
1253       CALL diffusion_s( pt, surf_def_h(0)%shf, surf_def_h(1)%shf, surf_def_h(2)%shf,              &
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,                  &
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 )
1260
1261!
1262!--    Consideration of heat sources within the plant canopy
1263       IF ( plant_canopy  .AND.  (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
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
1275       IF ( nudging )  CALL nudge( simulated_time, 'pt' )
1276
1277!
1278!--    If required compute influence of large-scale subsidence/ascent
1279       IF ( large_scale_subsidence  .AND.  .NOT. use_subsidence_tendencies )  THEN
1280          CALL subsidence( tend, pt, pt_init, 2 )
1281       ENDIF
1282
1283#if defined( __rrtmg )
1284!
1285!--    If required, add tendency due to radiative heating/cooling
1286       IF ( radiation  .AND.  simulated_time > skip_time_do_radiation )  THEN
1287            CALL radiation_tendency ( tend )
1288       ENDIF
1289#endif
1290
1291       CALL module_interface_actions( 'pt-tendency' )
1292
1293!
1294!--    Prognostic equation for potential temperature
1295       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1296       !$ACC PRESENT(pt, tend, tpt_m, wall_flags_total_0) &
1297       !$ACC PRESENT(pt_init, rdf_sc, ptdf_x, ptdf_y) &
1298       !$ACC PRESENT(tsc(3:5)) &
1299       !$ACC PRESENT(pt_p)
1300       DO  i = nxl, nxr
1301          DO  j = nys, nyn
1302!
1303!--          Following directive is required to vectorize on Intel19
1304             !DIR$ IVDEP
1305             DO  k = nzb+1, nzt
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 ) )
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
1318             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1319             !$ACC PRESENT(tend, tpt_m)
1320             DO  i = nxl, nxr
1321                DO  j = nys, nyn
1322                   DO  k = nzb+1, nzt
1323                      tpt_m(k,j,i) = tend(k,j,i)
1324                   ENDDO
1325                ENDDO
1326             ENDDO
1327          ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
1328             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1329             !$ACC PRESENT(tend, tpt_m)
1330             DO  i = nxl, nxr
1331                DO  j = nys, nyn
1332                   DO  k = nzb+1, nzt
1333                      tpt_m(k,j,i) = - 9.5625_wp * tend(k,j,i) + 5.3125_wp * tpt_m(k,j,i)
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!
1345!-- If required, compute prognostic equation for total water content
1346    IF ( humidity )  THEN
1347
1348       CALL cpu_log( log_point(29), 'q-equation', 'start' )
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
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,                             &
1375                                 bc_dirichlet_s  .OR.  bc_radiation_s )
1376             ELSE
1377                CALL advec_s_pw( q )
1378             ENDIF
1379          ELSE
1380             CALL advec_s_up( q )
1381          ENDIF
1382       ENDIF
1383
1384       CALL diffusion_s( q, surf_def_h(0)%qsws, surf_def_h(1)%qsws, surf_def_h(2)%qsws,            &
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,               &
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 )
1391
1392!
1393!--    Sink or source of humidity due to canopy elements
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
1404       IF ( nudging )  CALL nudge( simulated_time, 'q' )
1405
1406!
1407!--    If required compute influence of large-scale subsidence/ascent
1408       IF ( large_scale_subsidence  .AND.  .NOT. use_subsidence_tendencies )  THEN
1409         CALL subsidence( tend, q, q_init, 3 )
1410       ENDIF
1411
1412       CALL module_interface_actions( 'q-tendency' )
1413
1414!
1415!--    Prognostic equation for total water content
1416       DO  i = nxl, nxr
1417          DO  j = nys, nyn
1418!
1419!--          Following directive is required to vectorize on Intel19
1420             !DIR$ IVDEP
1421             DO  k = nzb+1, nzt
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 ) )
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
1437                   DO  k = nzb+1, nzt
1438                      tq_m(k,j,i) = tend(k,j,i)
1439                   ENDDO
1440                ENDDO
1441             ENDDO
1442          ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
1443             DO  i = nxl, nxr
1444                DO  j = nys, nyn
1445                   DO  k = nzb+1, nzt
1446                      tq_m(k,j,i) = - 9.5625_wp * tend(k,j,i) + 5.3125_wp * tq_m(k,j,i)
1447                   ENDDO
1448                ENDDO
1449             ENDDO
1450          ENDIF
1451       ENDIF
1452
1453       CALL cpu_log( log_point(29), 'q-equation', 'stop' )
1454
1455    ENDIF
1456!
1457!-- If required, compute prognostic equation for scalar
1458    IF ( passive_scalar )  THEN
1459
1460       CALL cpu_log( log_point(66), 's-equation', 'start' )
1461
1462!
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
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,                             &
1487                                 bc_dirichlet_s  .OR.  bc_radiation_s )
1488             ELSE
1489                CALL advec_s_pw( s )
1490             ENDIF
1491          ELSE
1492             CALL advec_s_up( s )
1493          ENDIF
1494       ENDIF
1495
1496       CALL diffusion_s( s, surf_def_h(0)%ssws, surf_def_h(1)%ssws, surf_def_h(2)%ssws,            &
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,               &
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 )
1503
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.
1516!        IF ( nudging )  CALL nudge( simulated_time, 'q' )
1517
1518!
1519!--    If required compute influence of large-scale subsidence/ascent.
1520!--    Not implemented for scalars so far.
1521       IF ( large_scale_subsidence  .AND.  .NOT. use_subsidence_tendencies  .AND.                  &
1522            .NOT. large_scale_forcing )  THEN
1523         CALL subsidence( tend, s, s_init, 3 )
1524       ENDIF
1525
1526       CALL module_interface_actions( 's-tendency' )
1527
1528!
1529!--    Prognostic equation for total water content
1530       DO  i = nxl, nxr
1531          DO  j = nys, nyn
1532!
1533!--          Following directive is required to vectorize on Intel19
1534             !DIR$ IVDEP
1535             DO  k = nzb+1, nzt
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 ) )
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
1551                   DO  k = nzb+1, nzt
1552                      ts_m(k,j,i) = tend(k,j,i)
1553                   ENDDO
1554                ENDDO
1555             ENDDO
1556          ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
1557             DO  i = nxl, nxr
1558                DO  j = nys, nyn
1559                   DO  k = nzb+1, nzt
1560                      ts_m(k,j,i) = - 9.5625_wp * tend(k,j,i) + 5.3125_wp * ts_m(k,j,i)
1561                   ENDDO
1562                ENDDO
1563             ENDDO
1564          ENDIF
1565       ENDIF
1566
1567       CALL cpu_log( log_point(66), 's-equation', 'stop' )
1568
1569    ENDIF
1570!
1571!-- Calculate prognostic equations for all other modules
1572    CALL module_interface_prognostic_equations()
1573
1574    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_vector', 'end' )
1575
1576 END SUBROUTINE prognostic_equations_vector
1577
1578
1579 END MODULE prognostic_equations_mod
Note: See TracBrowser for help on using the repository browser.