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

Last change on this file since 4360 was 4360, checked in by suehring, 16 months ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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