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

Last change on this file since 4652 was 4649, checked in by raasch, 4 years ago

files re-formatted to follow the PALM coding standard

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