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

Last change on this file since 4180 was 4180, checked in by scharf, 2 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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