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

Last change on this file since 4811 was 4731, checked in by schwenkel, 4 years ago

Move exchange_horiz from time_integration to modules

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