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

Last change on this file since 4598 was 4370, checked in by raasch, 16 months ago

bugfixes for previous commit: unused variables removed, Temperton-fft usage on GPU, openacc porting of vector version of Obukhov length calculation, collective read switched off on NEC to avoid hanging; some vector directives added in prognostic equations to force vectorization on Intel19 compiler, configuration files for NEC Aurora added

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