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

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