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

Last change on this file since 4048 was 4048, checked in by knoop, 2 years ago

Moved turbulence_closure_mod calls into module_interface

  • Property svn:keywords set to Id
File size: 63.7 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 4048 2019-06-21 21:00:21Z knoop $
27! Moved tcm_prognostic_equations to module_interface
28!
29! 3987 2019-05-22 09:52:13Z kanani
30! Introduce alternative switch for debug output during timestepping
31!
32! 3956 2019-05-07 12:32:52Z monakurppa
33! Removed salsa calls.
34!
35! 3931 2019-04-24 16:34:28Z schwenkel
36! Correct/complete module_interface introduction for chemistry model
37!
38! 3899 2019-04-16 14:05:27Z monakurppa
39! Corrections in the OpenMP version of salsa
40!
41! 3887 2019 -04-12 08:47:41Z schwenkel
42! Implicit Bugfix for chemistry model, loop for non_transport_physics over
43! ghost points is avoided. Instead introducing module_interface_exchange_horiz.
44!
45! 3885 2019-04-11 11:29:34Z kanani
46! Changes related to global restructuring of location messages and introduction
47! of additional debug messages
48!
49! 3881 2019-04-10 09:31:22Z suehring
50! Bugfix in OpenMP directive
51!
52! 3880 2019-04-08 21:43:02Z knoop
53! Moved wtm_tendencies to module_interface_actions
54!
55! 3874 2019-04-08 16:53:48Z knoop
56! Added non_transport_physics module interfaces and moved bcm code into it
57!
58! 3872 2019-04-08 15:03:06Z knoop
59! Moving prognostic equations of bcm into bulk_cloud_model_mod
60!
61! 3864 2019-04-05 09:01:56Z monakurppa
62! Modifications made for salsa:
63! - salsa_prognostic_equations moved to salsa_mod (and the call to
64!   module_interface_mod)
65! - Renamed nbins --> nbins_aerosol, ncc_tot --> ncomponents_mass and
66!   ngast --> ngases_salsa and loop indices b, c and sg to ib, ic and ig
67!
68! 3840 2019-03-29 10:35:52Z knoop
69! added USE chem_gasphase_mod for nspec, nspec and spc_names
70!
71! 3820 2019-03-27 11:53:41Z forkel
72! renamed do_depo to deposition_dry (ecc)
73!
74! 3797 2019-03-15 11:15:38Z forkel
75! Call chem_integegrate in OpenMP loop   (ketelsen)
76!
77!
78! 3771 2019-02-28 12:19:33Z raasch
79! preprocessor directivs fro rrtmg added
80!
81! 3761 2019-02-25 15:31:42Z raasch
82! unused variable removed
83!
84! 3719 2019-02-06 13:10:18Z kanani
85! Cleaned up chemistry cpu measurements
86!
87! 3684 2019-01-20 20:20:58Z knoop
88! OpenACC port for SPEC
89!
90! 3589 2018-11-30 15:09:51Z suehring
91! Move the control parameter "salsa" from salsa_mod to control_parameters
92! (M. Kurppa)
93!
94! 3582 2018-11-29 19:16:36Z suehring
95! Implementation of a new aerosol module salsa.
96! Remove cpu-logs from i,j loop in cache version.
97!
98! 3458 2018-10-30 14:51:23Z kanani
99! remove duplicate USE chem_modules
100! from chemistry branch r3443, banzhafs, basit:
101! chem_depo call introduced
102! code added for decycling chemistry
103!
104! 3386 2018-10-19 16:28:22Z gronemeier
105! Renamed tcm_prognostic to tcm_prognostic_equations
106!
107! 3355 2018-10-16 14:03:34Z knoop
108! (from branch resler)
109! Fix for chemistry call
110!
111! 3302 2018-10-03 02:39:40Z raasch
112! Stokes drift + wave breaking term added
113!
114! 3298 2018-10-02 12:21:11Z kanani
115! Code added for decycling chemistry (basit)
116!
117! 3294 2018-10-01 02:37:10Z raasch
118! changes concerning modularization of ocean option
119!
120! 3274 2018-09-24 15:42:55Z knoop
121! Modularization of all bulk cloud physics code components
122!
123! 3241 2018-09-12 15:02:00Z raasch
124! omp_get_thread_num now declared in omp directive
125!
126! 3183 2018-07-27 14:25:55Z suehring
127! Remove unused variables from USE statements
128!
129! 3182 2018-07-27 13:36:03Z suehring
130! Revise recent bugfix for nesting
131!
132! 3021 2018-05-16 08:14:20Z maronga
133! Bugfix in IF clause for nesting
134!
135! 3014 2018-05-09 08:42:38Z maronga
136! Fixed a bug in the IF condition to call pcm_tendency in case of
137! potential temperature
138!
139! 2815 2018-02-19 11:29:57Z kanani
140! Rename chem_tendency to chem_prognostic_equations,
141! implement vector version for air chemistry
142!
143! 2766 2018-01-22 17:17:47Z kanani
144! Removed preprocessor directive __chem
145!
146! 2746 2018-01-15 12:06:04Z suehring
147! Move flag plant canopy to modules
148!
149! 2719 2018-01-02 09:02:06Z maronga
150! Bugfix for last change.
151!
152! 2718 2018-01-02 08:49:38Z maronga
153! Corrected "Former revisions" section
154!
155! 2696 2017-12-14 17:12:51Z kanani
156! - Change in file header (GPL part)
157! - Moved TKE equation to tcm_prognostic (TG)
158! - Added switch for chemical reactions (RF, FK)
159! - Implementation of chemistry module (RF, BK, FK)
160!
161! 2563 2017-10-19 15:36:10Z Giersch
162! Variable wind_turbine moved to module control_parameters
163!
164! 2320 2017-07-21 12:47:43Z suehring
165! Modularize large-scale forcing and nudging
166!
167! 2292 2017-06-20 09:51:42Z schwenkel
168! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
169! includes two more prognostic equations for cloud drop concentration (nc) 
170! and cloud water content (qc).
171!
172! 2261 2017-06-08 14:25:57Z raasch
173! bugfix for r2232: openmp directives removed
174!
175! 2233 2017-05-30 18:08:54Z suehring
176!
177! 2232 2017-05-30 17:47:52Z suehring
178! Adjutst to new surface-type structure. Remove call for usm_wall_heat_flux,
179! which is realized directly in diffusion_s now.
180!
181! 2192 2017-03-22 04:14:10Z raasch
182! Bugfix for misplaced and missing openMP directives from r2155
183!
184! 2155 2017-02-21 09:57:40Z hoffmann
185! Bugfix in the calculation of microphysical quantities on ghost points.
186!
187! 2118 2017-01-17 16:38:49Z raasch
188! OpenACC version of subroutine removed
189!
190! 2031 2016-10-21 15:11:58Z knoop
191! renamed variable rho to rho_ocean
192!
193! 2011 2016-09-19 17:29:57Z kanani
194! Flag urban_surface is now defined in module control_parameters.
195!
196! 2007 2016-08-24 15:47:17Z kanani
197! Added pt tendency calculation based on energy balance at urban surfaces
198! (new urban surface model)
199!
200! 2000 2016-08-20 18:09:15Z knoop
201! Forced header and separation lines into 80 columns
202!
203! 1976 2016-07-27 13:28:04Z maronga
204! Simplied calls to radiation model
205!
206! 1960 2016-07-12 16:34:24Z suehring
207! Separate humidity and passive scalar
208!
209! 1914 2016-05-26 14:44:07Z witha
210! Added calls for wind turbine model
211!
212! 1873 2016-04-18 14:50:06Z maronga
213! Module renamed (removed _mod)
214!
215! 1850 2016-04-08 13:29:27Z maronga
216! Module renamed
217!
218! 1826 2016-04-07 12:01:39Z maronga
219! Renamed canopy model calls.
220!
221! 1822 2016-04-07 07:49:42Z hoffmann
222! Kessler microphysics scheme moved to microphysics.
223!
224! 1757 2016-02-22 15:49:32Z maronga
225!
226! 1691 2015-10-26 16:17:44Z maronga
227! Added optional model spin-up without radiation / land surface model calls.
228! Formatting corrections.
229!
230! 1682 2015-10-07 23:56:08Z knoop
231! Code annotations made doxygen readable
232!
233! 1585 2015-04-30 07:05:52Z maronga
234! Added call for temperature tendency calculation due to radiative flux divergence
235!
236! 1517 2015-01-07 19:12:25Z hoffmann
237! advec_s_bc_mod addded, since advec_s_bc is now a module
238!
239! 1484 2014-10-21 10:53:05Z kanani
240! Changes due to new module structure of the plant canopy model:
241! parameters cthf and plant_canopy moved to module plant_canopy_model_mod.
242! Removed double-listing of use_upstream_for_tke in ONLY-list of module
243! control_parameters
244!
245! 1409 2014-05-23 12:11:32Z suehring
246! Bugfix: i_omp_start changed for advec_u_ws at left inflow and outflow boundary.
247! This ensures that left-hand side fluxes are also calculated for nxl in that
248! case, even though the solution at nxl is overwritten in boundary_conds()
249!
250! 1398 2014-05-07 11:15:00Z heinze
251! Rayleigh-damping for horizontal velocity components changed: instead of damping
252! against ug and vg, damping against u_init and v_init is used to allow for a
253! homogenized treatment in case of nudging
254!
255! 1380 2014-04-28 12:40:45Z heinze
256! Change order of calls for scalar prognostic quantities:
257! ls_advec -> nudging -> subsidence since initial profiles
258!
259! 1374 2014-04-25 12:55:07Z raasch
260! missing variables added to ONLY lists
261!
262! 1365 2014-04-22 15:03:56Z boeske
263! Calls of ls_advec for large scale advection added,
264! subroutine subsidence is only called if use_subsidence_tendencies = .F.,
265! new argument ls_index added to the calls of subsidence
266! +ls_index
267!
268! 1361 2014-04-16 15:17:48Z hoffmann
269! Two-moment microphysics moved to the start of prognostic equations. This makes
270! the 3d arrays for tend_q, tend_qr, tend_pt and tend_pt redundant.
271! Additionally, it is allowed to call the microphysics just once during the time
272! step (not at each sub-time step).
273!
274! Two-moment cloud physics added for vector and accelerator optimization.
275!
276! 1353 2014-04-08 15:21:23Z heinze
277! REAL constants provided with KIND-attribute
278!
279! 1337 2014-03-25 15:11:48Z heinze
280! Bugfix: REAL constants provided with KIND-attribute
281!
282! 1332 2014-03-25 11:59:43Z suehring
283! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke
284!
285! 1330 2014-03-24 17:29:32Z suehring
286! In case of SGS-particle velocity advection of TKE is also allowed with
287! dissipative 5th-order scheme.
288!
289! 1320 2014-03-20 08:40:49Z raasch
290! ONLY-attribute added to USE-statements,
291! kind-parameters added to all INTEGER and REAL declaration statements,
292! kinds are defined in new module kinds,
293! old module precision_kind is removed,
294! revision history before 2012 removed,
295! comment fields (!:) to be used for variable explanations added to
296! all variable declaration statements
297!
298! 1318 2014-03-17 13:35:16Z raasch
299! module interfaces removed
300!
301! 1257 2013-11-08 15:18:40Z raasch
302! openacc loop vector clauses removed, independent clauses added
303!
304! 1246 2013-11-01 08:59:45Z heinze
305! enable nudging also for accelerator version
306!
307! 1241 2013-10-30 11:36:58Z heinze
308! usage of nudging enabled (so far not implemented for accelerator version)
309!
310! 1179 2013-06-14 05:57:58Z raasch
311! two arguments removed from routine buoyancy, ref_state updated on device
312!
313! 1128 2013-04-12 06:19:32Z raasch
314! those parts requiring global communication moved to time_integration,
315! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
316! j_north
317!
318! 1115 2013-03-26 18:16:16Z hoffmann
319! optimized cloud physics: calculation of microphysical tendencies transfered
320! to microphysics.f90; qr and nr are only calculated if precipitation is required
321!
322! 1111 2013-03-08 23:54:10Z raasch
323! update directives for prognostic quantities removed
324!
325! 1106 2013-03-04 05:31:38Z raasch
326! small changes in code formatting
327!
328! 1092 2013-02-02 11:24:22Z raasch
329! unused variables removed
330!
331! 1053 2012-11-13 17:11:03Z hoffmann
332! implementation of two new prognostic equations for rain drop concentration (nr)
333! and rain water content (qr)
334!
335! currently, only available for cache loop optimization
336!
337! 1036 2012-10-22 13:43:42Z raasch
338! code put under GPL (PALM 3.9)
339!
340! 1019 2012-09-28 06:46:45Z raasch
341! non-optimized version of prognostic_equations removed
342!
343! 1015 2012-09-27 09:23:24Z raasch
344! new branch prognostic_equations_acc
345! OpenACC statements added + code changes required for GPU optimization
346!
347! 1001 2012-09-13 14:08:46Z raasch
348! all actions concerning leapfrog- and upstream-spline-scheme removed
349!
350! 978 2012-08-09 08:28:32Z fricke
351! km_damp_x and km_damp_y removed in calls of diffusion_u and diffusion_v
352! add ptdf_x, ptdf_y for damping the potential temperature at the inflow
353! boundary in case of non-cyclic lateral boundaries
354! Bugfix: first thread index changes for WS-scheme at the inflow
355!
356! 940 2012-07-09 14:31:00Z raasch
357! temperature equation can be switched off
358!
359! Revision 1.1  2000/04/13 14:56:27  schroeter
360! Initial revision
361!
362!
363! Description:
364! ------------
365!> Solving the prognostic equations.
366!------------------------------------------------------------------------------!
367 MODULE prognostic_equations_mod
368
369    USE advec_s_bc_mod,                                                        &
370        ONLY:  advec_s_bc
371
372    USE advec_s_pw_mod,                                                        &
373        ONLY:  advec_s_pw
374
375    USE advec_s_up_mod,                                                        &
376        ONLY:  advec_s_up
377
378    USE advec_u_pw_mod,                                                        &
379        ONLY:  advec_u_pw
380
381    USE advec_u_up_mod,                                                        &
382        ONLY:  advec_u_up
383
384    USE advec_v_pw_mod,                                                        &
385        ONLY:  advec_v_pw
386
387    USE advec_v_up_mod,                                                        &
388        ONLY:  advec_v_up
389
390    USE advec_w_pw_mod,                                                        &
391        ONLY:  advec_w_pw
392
393    USE advec_w_up_mod,                                                        &
394        ONLY:  advec_w_up
395
396    USE advec_ws,                                                              &
397        ONLY:  advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws
398
399    USE arrays_3d,                                                             &
400        ONLY:  diss_l_e, diss_l_pt, diss_l_q,                                  &
401               diss_l_s, diss_l_sa, diss_s_e,                                  &
402               diss_s_pt, diss_s_q, diss_s_s, diss_s_sa,                       &
403               e, e_p, flux_s_e, flux_s_pt, flux_s_q,                          &
404               flux_s_s, flux_s_sa, flux_l_e,                                  &
405               flux_l_pt, flux_l_q, flux_l_s,                                  &
406               flux_l_sa, pt, ptdf_x, ptdf_y, pt_init,                         &
407               pt_p, prho, q, q_init, q_p, rdf, rdf_sc,                        &
408               ref_state, rho_ocean, s, s_init, s_p, tend, te_m,               &
409               tpt_m, tq_m, ts_m, tu_m, tv_m, tw_m, u,                         &
410               ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
411
412    USE buoyancy_mod,                                                          &
413        ONLY:  buoyancy
414
415    USE control_parameters,                                                    &
416        ONLY:  constant_diffusion,                                             &
417               debug_output_timestep,                                          &
418               dp_external, dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d,    &
419               humidity, intermediate_timestep_count,                          &
420               intermediate_timestep_count_max, large_scale_forcing,           &
421               large_scale_subsidence, neutral, nudging,                       &
422               ocean_mode, passive_scalar, plant_canopy, pt_reference,         &
423               scalar_advec, scalar_advec, simulated_time, sloping_surface,    &
424               timestep_scheme, tsc, use_subsidence_tendencies,                &
425               use_upstream_for_tke, wind_turbine, ws_scheme_mom,              & 
426               ws_scheme_sca, urban_surface, land_surface,                     &
427               time_since_reference_point, salsa
428
429    USE coriolis_mod,                                                          &
430        ONLY:  coriolis
431
432    USE cpulog,                                                                &
433        ONLY:  cpu_log, log_point, log_point_s
434
435    USE diffusion_s_mod,                                                       &
436        ONLY:  diffusion_s
437
438    USE diffusion_u_mod,                                                       &
439        ONLY:  diffusion_u
440
441    USE diffusion_v_mod,                                                       &
442        ONLY:  diffusion_v
443
444    USE diffusion_w_mod,                                                       &
445        ONLY:  diffusion_w
446
447    USE indices,                                                               &
448        ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
449               nzb, nzt, wall_flags_0
450
451    USE kinds
452
453    USE lsf_nudging_mod,                                                       &
454        ONLY:  ls_advec, nudge
455
456    USE module_interface,                                                      &
457        ONLY:  module_interface_actions, &
458               module_interface_non_advective_processes, &
459               module_interface_exchange_horiz, &
460               module_interface_prognostic_equations
461
462    USE ocean_mod,                                                             &
463        ONLY:  stokes_drift_terms, stokes_force,   &
464               wave_breaking, wave_breaking_term
465
466    USE plant_canopy_model_mod,                                                &
467        ONLY:  cthf, pcm_tendency
468
469#if defined( __rrtmg )
470    USE radiation_model_mod,                                                   &
471        ONLY:  radiation, radiation_tendency,                                  &
472               skip_time_do_radiation
473#endif
474
475    USE statistics,                                                            &
476        ONLY:  hom
477
478    USE subsidence_mod,                                                        &
479        ONLY:  subsidence
480
481    USE surface_mod,                                                           &
482        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
483                surf_usm_v
484
485    IMPLICIT NONE
486
487    PRIVATE
488    PUBLIC prognostic_equations_cache, prognostic_equations_vector
489
490    INTERFACE prognostic_equations_cache
491       MODULE PROCEDURE prognostic_equations_cache
492    END INTERFACE prognostic_equations_cache
493
494    INTERFACE prognostic_equations_vector
495       MODULE PROCEDURE prognostic_equations_vector
496    END INTERFACE prognostic_equations_vector
497
498
499 CONTAINS
500
501
502!------------------------------------------------------------------------------!
503! Description:
504! ------------
505!> Version with one optimized loop over all equations. It is only allowed to
506!> be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
507!>
508!> Here the calls of most subroutines are embedded in two DO loops over i and j,
509!> so communication between CPUs is not allowed (does not make sense) within
510!> these loops.
511!>
512!> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
513!------------------------------------------------------------------------------!
514
515 SUBROUTINE prognostic_equations_cache
516
517
518    INTEGER(iwp) ::  i                   !<
519    INTEGER(iwp) ::  i_omp_start         !<
520    INTEGER(iwp) ::  j                   !<
521    INTEGER(iwp) ::  k                   !<
522!$  INTEGER(iwp) ::  omp_get_thread_num  !<
523    INTEGER(iwp) ::  tn = 0              !<
524
525    LOGICAL      ::  loop_start          !<
526
527
528    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_cache', 'start' )
529!
530!-- Time measurement can only be performed for the whole set of equations
531    CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
532
533    !$OMP PARALLEL PRIVATE (i,j)
534    !$OMP DO
535    DO  i = nxl, nxr
536       DO  j = nys, nyn
537!
538!--       Calculate non advective processes for all other modules
539          CALL module_interface_non_advective_processes( i, j )
540       ENDDO
541    ENDDO
542!
543!-- Module Inferface for exchange horiz after non_advective_processes but before
544!-- advection. Therefore, non_advective_processes must not run for ghost points.
545    !$OMP END PARALLEL
546    CALL module_interface_exchange_horiz()
547!
548!-- Loop over all prognostic equations
549    !$OMP PARALLEL PRIVATE (i,i_omp_start,j,k,loop_start,tn)
550
551    !$  tn = omp_get_thread_num()
552    loop_start = .TRUE.
553
554    !$OMP DO
555    DO  i = nxl, nxr
556
557!
558!--    Store the first loop index. It differs for each thread and is required
559!--    later in advec_ws
560       IF ( loop_start )  THEN
561          loop_start  = .FALSE.
562          i_omp_start = i
563       ENDIF
564
565       DO  j = nys, nyn
566!
567!--       Tendency terms for u-velocity component. Please note, in case of
568!--       non-cyclic boundary conditions the grid point i=0 is excluded from
569!--       the prognostic equations for the u-component.
570          IF ( i >= nxlu )  THEN
571
572             tend(:,j,i) = 0.0_wp
573             IF ( timestep_scheme(1:5) == 'runge' )  THEN
574                IF ( ws_scheme_mom )  THEN
575                   CALL advec_u_ws( i, j, i_omp_start, tn )
576                ELSE
577                   CALL advec_u_pw( i, j )
578                ENDIF
579             ELSE
580                CALL advec_u_up( i, j )
581             ENDIF
582             CALL diffusion_u( i, j )
583             CALL coriolis( i, j, 1 )
584             IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
585                CALL buoyancy( i, j, pt, 1 )
586             ENDIF
587
588!
589!--          Drag by plant canopy
590             IF ( plant_canopy )  CALL pcm_tendency( i, j, 1 )
591
592!
593!--          External pressure gradient
594             IF ( dp_external )  THEN
595                DO  k = dp_level_ind_b+1, nzt
596                   tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
597                ENDDO
598             ENDIF
599
600!
601!--          Nudging
602             IF ( nudging )  CALL nudge( i, j, simulated_time, 'u' )
603
604!
605!--          Effect of Stokes drift (in ocean mode only)
606             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 1 )
607
608             CALL module_interface_actions( i, j, 'u-tendency' )
609!
610!--          Prognostic equation for u-velocity component
611             DO  k = nzb+1, nzt
612
613                u_p(k,j,i) = u(k,j,i) + ( dt_3d *                               &
614                                            ( tsc(2) * tend(k,j,i) +            &
615                                              tsc(3) * tu_m(k,j,i) )            &
616                                            - tsc(5) * rdf(k)                   &
617                                                     * ( u(k,j,i) - u_init(k) ) &
618                                        ) * MERGE( 1.0_wp, 0.0_wp,              &
619                                                   BTEST( wall_flags_0(k,j,i), 1 )&
620                                                 )
621             ENDDO
622
623!
624!--          Add turbulence generated by wave breaking (in ocean mode only)
625             IF ( wave_breaking  .AND.                                         &
626               intermediate_timestep_count == intermediate_timestep_count_max )&
627             THEN
628                CALL wave_breaking_term( i, j, 1 )
629             ENDIF
630
631!
632!--          Calculate tendencies for the next Runge-Kutta step
633             IF ( timestep_scheme(1:5) == 'runge' )  THEN
634                IF ( intermediate_timestep_count == 1 )  THEN
635                   DO  k = nzb+1, nzt
636                      tu_m(k,j,i) = tend(k,j,i)
637                   ENDDO
638                ELSEIF ( intermediate_timestep_count < &
639                         intermediate_timestep_count_max )  THEN
640                   DO  k = nzb+1, nzt
641                      tu_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                &
642                                     + 5.3125_wp * tu_m(k,j,i)
643                   ENDDO
644                ENDIF
645             ENDIF
646
647          ENDIF
648!
649!--       Tendency terms for v-velocity component. Please note, in case of
650!--       non-cyclic boundary conditions the grid point j=0 is excluded from
651!--       the prognostic equations for the v-component. !--       
652          IF ( j >= nysv )  THEN
653
654             tend(:,j,i) = 0.0_wp
655             IF ( timestep_scheme(1:5) == 'runge' )  THEN
656                IF ( ws_scheme_mom )  THEN
657                    CALL advec_v_ws( i, j, i_omp_start, tn )
658                ELSE
659                    CALL advec_v_pw( i, j )
660                ENDIF
661             ELSE
662                CALL advec_v_up( i, j )
663             ENDIF
664             CALL diffusion_v( i, j )
665             CALL coriolis( i, j, 2 )
666
667!
668!--          Drag by plant canopy
669             IF ( plant_canopy )  CALL pcm_tendency( i, j, 2 )
670
671!
672!--          External pressure gradient
673             IF ( dp_external )  THEN
674                DO  k = dp_level_ind_b+1, nzt
675                   tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
676                ENDDO
677             ENDIF
678
679!
680!--          Nudging
681             IF ( nudging )  CALL nudge( i, j, simulated_time, 'v' )
682
683!
684!--          Effect of Stokes drift (in ocean mode only)
685             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 2 )
686
687             CALL module_interface_actions( i, j, 'v-tendency' )
688!
689!--          Prognostic equation for v-velocity component
690             DO  k = nzb+1, nzt
691                v_p(k,j,i) = v(k,j,i) + ( dt_3d *                              &
692                                            ( tsc(2) * tend(k,j,i) +           &
693                                              tsc(3) * tv_m(k,j,i) )           &
694                                            - tsc(5) * rdf(k)                  &
695                                                     * ( v(k,j,i) - v_init(k) )&
696                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
697                                                   BTEST( wall_flags_0(k,j,i), 2 )&
698                                                 )
699             ENDDO
700
701!
702!--          Add turbulence generated by wave breaking (in ocean mode only)
703             IF ( wave_breaking  .AND.                                         &
704               intermediate_timestep_count == intermediate_timestep_count_max )&
705             THEN
706                CALL wave_breaking_term( i, j, 2 )
707             ENDIF
708
709!
710!--          Calculate tendencies for the next Runge-Kutta step
711             IF ( timestep_scheme(1:5) == 'runge' )  THEN
712                IF ( intermediate_timestep_count == 1 )  THEN
713                   DO  k = nzb+1, nzt
714                      tv_m(k,j,i) = tend(k,j,i)
715                   ENDDO
716                ELSEIF ( intermediate_timestep_count < &
717                         intermediate_timestep_count_max )  THEN
718                   DO  k = nzb+1, nzt
719                      tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
720                                     + 5.3125_wp * tv_m(k,j,i)
721                   ENDDO
722                ENDIF
723             ENDIF
724
725          ENDIF
726
727!
728!--       Tendency terms for w-velocity component
729          tend(:,j,i) = 0.0_wp
730          IF ( timestep_scheme(1:5) == 'runge' )  THEN
731             IF ( ws_scheme_mom )  THEN
732                CALL advec_w_ws( i, j, i_omp_start, tn )
733             ELSE
734                CALL advec_w_pw( i, j )
735             END IF
736          ELSE
737             CALL advec_w_up( i, j )
738          ENDIF
739          CALL diffusion_w( i, j )
740          CALL coriolis( i, j, 3 )
741
742          IF ( .NOT. neutral )  THEN
743             IF ( ocean_mode )  THEN
744                CALL buoyancy( i, j, rho_ocean, 3 )
745             ELSE
746                IF ( .NOT. humidity )  THEN
747                   CALL buoyancy( i, j, pt, 3 )
748                ELSE
749                   CALL buoyancy( i, j, vpt, 3 )
750                ENDIF
751             ENDIF
752          ENDIF
753
754!
755!--       Drag by plant canopy
756          IF ( plant_canopy )  CALL pcm_tendency( i, j, 3 )
757
758!
759!--       Effect of Stokes drift (in ocean mode only)
760          IF ( stokes_force )  CALL stokes_drift_terms( i, j, 3 )
761
762          CALL module_interface_actions( i, j, 'w-tendency' )
763!
764!--       Prognostic equation for w-velocity component
765          DO  k = nzb+1, nzt-1
766             w_p(k,j,i) = w(k,j,i) + ( dt_3d *                                 &
767                                             ( tsc(2) * tend(k,j,i) +          &
768                                               tsc(3) * tw_m(k,j,i) )          &
769                                             - tsc(5) * rdf(k) * w(k,j,i)      &
770                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
771                                                BTEST( wall_flags_0(k,j,i), 3 )&
772                                              )
773          ENDDO
774
775!
776!--       Calculate tendencies for the next Runge-Kutta step
777          IF ( timestep_scheme(1:5) == 'runge' )  THEN
778             IF ( intermediate_timestep_count == 1 )  THEN
779                DO  k = nzb+1, nzt-1
780                   tw_m(k,j,i) = tend(k,j,i)
781                ENDDO
782             ELSEIF ( intermediate_timestep_count < &
783                      intermediate_timestep_count_max )  THEN
784                DO  k = nzb+1, nzt-1
785                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
786                                  + 5.3125_wp * tw_m(k,j,i)
787                ENDDO
788             ENDIF
789          ENDIF
790
791!
792!--       If required, compute prognostic equation for potential temperature
793          IF ( .NOT. neutral )  THEN
794!
795!--          Tendency terms for potential temperature
796             tend(:,j,i) = 0.0_wp
797             IF ( timestep_scheme(1:5) == 'runge' )  THEN
798                   IF ( ws_scheme_sca )  THEN
799                      CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, diss_s_pt,   &
800                                       flux_l_pt, diss_l_pt, i_omp_start, tn )
801                   ELSE
802                      CALL advec_s_pw( i, j, pt )
803                   ENDIF
804             ELSE
805                CALL advec_s_up( i, j, pt )
806             ENDIF
807             CALL diffusion_s( i, j, pt,                                       &
808                               surf_def_h(0)%shf, surf_def_h(1)%shf,           &
809                               surf_def_h(2)%shf,                              &
810                               surf_lsm_h%shf,    surf_usm_h%shf,              &
811                               surf_def_v(0)%shf, surf_def_v(1)%shf,           &
812                               surf_def_v(2)%shf, surf_def_v(3)%shf,           &
813                               surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,           &
814                               surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,           &
815                               surf_usm_v(0)%shf, surf_usm_v(1)%shf,           &
816                               surf_usm_v(2)%shf, surf_usm_v(3)%shf )
817
818!
819!--          Consideration of heat sources within the plant canopy
820             IF ( plant_canopy  .AND.                                          &
821                (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
822                CALL pcm_tendency( i, j, 4 )
823             ENDIF
824
825!
826!--          Large scale advection
827             IF ( large_scale_forcing )  THEN
828                CALL ls_advec( i, j, simulated_time, 'pt' )
829             ENDIF
830
831!
832!--          Nudging
833             IF ( nudging )  CALL nudge( i, j, simulated_time, 'pt' )
834
835!
836!--          If required, compute effect of large-scale subsidence/ascent
837             IF ( large_scale_subsidence  .AND.                                &
838                  .NOT. use_subsidence_tendencies )  THEN
839                CALL subsidence( i, j, tend, pt, pt_init, 2 )
840             ENDIF
841
842#if defined( __rrtmg )
843!
844!--          If required, add tendency due to radiative heating/cooling
845             IF ( radiation  .AND.                                             &
846                  simulated_time > skip_time_do_radiation )  THEN
847                CALL radiation_tendency ( i, j, tend )
848             ENDIF
849#endif
850
851             CALL module_interface_actions( i, j, 'pt-tendency' )
852!
853!--          Prognostic equation for potential temperature
854             DO  k = nzb+1, nzt
855                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d *                            &
856                                                  ( tsc(2) * tend(k,j,i) +     &
857                                                    tsc(3) * tpt_m(k,j,i) )    &
858                                                  - tsc(5)                     &
859                                                  * ( pt(k,j,i) - pt_init(k) ) &
860                                                  * ( rdf_sc(k) + ptdf_x(i)    &
861                                                                + ptdf_y(j) )  &
862                                          )                                    &
863                                       * MERGE( 1.0_wp, 0.0_wp,                &
864                                                BTEST( wall_flags_0(k,j,i), 0 )&
865                                              )
866             ENDDO
867
868!
869!--          Calculate tendencies for the next Runge-Kutta step
870             IF ( timestep_scheme(1:5) == 'runge' )  THEN
871                IF ( intermediate_timestep_count == 1 )  THEN
872                   DO  k = nzb+1, nzt
873                      tpt_m(k,j,i) = tend(k,j,i)
874                   ENDDO
875                ELSEIF ( intermediate_timestep_count < &
876                         intermediate_timestep_count_max )  THEN
877                   DO  k = nzb+1, nzt
878                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
879                                        5.3125_wp * tpt_m(k,j,i)
880                   ENDDO
881                ENDIF
882             ENDIF
883
884          ENDIF
885
886!
887!--       If required, compute prognostic equation for total water content
888          IF ( humidity )  THEN
889
890!
891!--          Tendency-terms for total water content / scalar
892             tend(:,j,i) = 0.0_wp
893             IF ( timestep_scheme(1:5) == 'runge' ) &
894             THEN
895                IF ( ws_scheme_sca )  THEN
896                   CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
897                                diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
898                ELSE
899                   CALL advec_s_pw( i, j, q )
900                ENDIF
901             ELSE
902                CALL advec_s_up( i, j, q )
903             ENDIF
904             CALL diffusion_s( i, j, q,                                        &
905                               surf_def_h(0)%qsws, surf_def_h(1)%qsws,         &
906                               surf_def_h(2)%qsws,                             &
907                               surf_lsm_h%qsws,    surf_usm_h%qsws,            &
908                               surf_def_v(0)%qsws, surf_def_v(1)%qsws,         &
909                               surf_def_v(2)%qsws, surf_def_v(3)%qsws,         &
910                               surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,         &
911                               surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,         &
912                               surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,         &
913                               surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
914
915!
916!--          Sink or source of humidity due to canopy elements
917             IF ( plant_canopy )  CALL pcm_tendency( i, j, 5 )
918
919!
920!--          Large scale advection
921             IF ( large_scale_forcing )  THEN
922                CALL ls_advec( i, j, simulated_time, 'q' )
923             ENDIF
924
925!
926!--          Nudging
927             IF ( nudging )  CALL nudge( i, j, simulated_time, 'q' )
928
929!
930!--          If required compute influence of large-scale subsidence/ascent
931             IF ( large_scale_subsidence  .AND.                                &
932                  .NOT. use_subsidence_tendencies )  THEN
933                CALL subsidence( i, j, tend, q, q_init, 3 )
934             ENDIF
935
936             CALL module_interface_actions( i, j, 'q-tendency' )
937
938!
939!--          Prognostic equation for total water content / scalar
940             DO  k = nzb+1, nzt
941                q_p(k,j,i) = q(k,j,i) + ( dt_3d *                              &
942                                                ( tsc(2) * tend(k,j,i) +       &
943                                                  tsc(3) * tq_m(k,j,i) )       &
944                                                - tsc(5) * rdf_sc(k) *         &
945                                                      ( q(k,j,i) - q_init(k) ) &
946                                        )                                      &
947                                       * MERGE( 1.0_wp, 0.0_wp,                &
948                                                BTEST( wall_flags_0(k,j,i), 0 )&
949                                              )               
950                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
951             ENDDO
952
953!
954!--          Calculate tendencies for the next Runge-Kutta step
955             IF ( timestep_scheme(1:5) == 'runge' )  THEN
956                IF ( intermediate_timestep_count == 1 )  THEN
957                   DO  k = nzb+1, nzt
958                      tq_m(k,j,i) = tend(k,j,i)
959                   ENDDO
960                ELSEIF ( intermediate_timestep_count < &
961                         intermediate_timestep_count_max )  THEN
962                   DO  k = nzb+1, nzt
963                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
964                                       5.3125_wp * tq_m(k,j,i)
965                   ENDDO
966                ENDIF
967             ENDIF
968
969          ENDIF
970
971!
972!--       If required, compute prognostic equation for scalar
973          IF ( passive_scalar )  THEN
974!
975!--          Tendency-terms for total water content / scalar
976             tend(:,j,i) = 0.0_wp
977             IF ( timestep_scheme(1:5) == 'runge' ) &
978             THEN
979                IF ( ws_scheme_sca )  THEN
980                   CALL advec_s_ws( i, j, s, 's', flux_s_s, &
981                                diss_s_s, flux_l_s, diss_l_s, i_omp_start, tn )
982                ELSE
983                   CALL advec_s_pw( i, j, s )
984                ENDIF
985             ELSE
986                CALL advec_s_up( i, j, s )
987             ENDIF
988             CALL diffusion_s( i, j, s,                                        &
989                               surf_def_h(0)%ssws, surf_def_h(1)%ssws,         &
990                               surf_def_h(2)%ssws,                             &
991                               surf_lsm_h%ssws,    surf_usm_h%ssws,            &
992                               surf_def_v(0)%ssws, surf_def_v(1)%ssws,         &
993                               surf_def_v(2)%ssws, surf_def_v(3)%ssws,         &
994                               surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,         &
995                               surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,         &
996                               surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,         &
997                               surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
998
999!
1000!--          Sink or source of scalar concentration due to canopy elements
1001             IF ( plant_canopy )  CALL pcm_tendency( i, j, 7 )
1002
1003!
1004!--          Large scale advection, still need to be extended for scalars
1005!              IF ( large_scale_forcing )  THEN
1006!                 CALL ls_advec( i, j, simulated_time, 's' )
1007!              ENDIF
1008
1009!
1010!--          Nudging, still need to be extended for scalars
1011!              IF ( nudging )  CALL nudge( i, j, simulated_time, 's' )
1012
1013!
1014!--          If required compute influence of large-scale subsidence/ascent.
1015!--          Note, the last argument is of no meaning in this case, as it is
1016!--          only used in conjunction with large_scale_forcing, which is to
1017!--          date not implemented for scalars.
1018             IF ( large_scale_subsidence  .AND.                                &
1019                  .NOT. use_subsidence_tendencies  .AND.                       &
1020                  .NOT. large_scale_forcing )  THEN
1021                CALL subsidence( i, j, tend, s, s_init, 3 )
1022             ENDIF
1023
1024             CALL module_interface_actions( i, j, 's-tendency' )
1025
1026!
1027!--          Prognostic equation for scalar
1028             DO  k = nzb+1, nzt
1029                s_p(k,j,i) = s(k,j,i) + (  dt_3d *                             &
1030                                            ( tsc(2) * tend(k,j,i) +           &
1031                                              tsc(3) * ts_m(k,j,i) )           &
1032                                            - tsc(5) * rdf_sc(k)               &
1033                                                     * ( s(k,j,i) - s_init(k) )&
1034                                        )                                      &
1035                                       * MERGE( 1.0_wp, 0.0_wp,                &
1036                                                BTEST( wall_flags_0(k,j,i), 0 )&
1037                                              )
1038                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
1039             ENDDO
1040
1041!
1042!--          Calculate tendencies for the next Runge-Kutta step
1043             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1044                IF ( intermediate_timestep_count == 1 )  THEN
1045                   DO  k = nzb+1, nzt
1046                      ts_m(k,j,i) = tend(k,j,i)
1047                   ENDDO
1048                ELSEIF ( intermediate_timestep_count < &
1049                         intermediate_timestep_count_max )  THEN
1050                   DO  k = nzb+1, nzt
1051                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
1052                                       5.3125_wp * ts_m(k,j,i)
1053                   ENDDO
1054                ENDIF
1055             ENDIF
1056
1057          ENDIF
1058!
1059!--       Calculate prognostic equations for all other modules
1060          CALL module_interface_prognostic_equations( i, j, i_omp_start, tn )
1061
1062       ENDDO  ! loop over j
1063    ENDDO  ! loop over i
1064    !$OMP END PARALLEL
1065
1066
1067    CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1068
1069    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_cache', 'end' )
1070
1071 END SUBROUTINE prognostic_equations_cache
1072
1073
1074!------------------------------------------------------------------------------!
1075! Description:
1076! ------------
1077!> Version for vector machines
1078!------------------------------------------------------------------------------!
1079
1080 SUBROUTINE prognostic_equations_vector
1081
1082
1083    INTEGER(iwp) ::  i     !<
1084    INTEGER(iwp) ::  j     !<
1085    INTEGER(iwp) ::  k     !<
1086
1087    REAL(wp)     ::  sbt  !<
1088
1089
1090    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_vector', 'start' )
1091!
1092!-- Calculate non advective processes for all other modules
1093    CALL module_interface_non_advective_processes
1094!
1095!-- Module Inferface for exchange horiz after non_advective_processes but before
1096!-- advection. Therefore, non_advective_processes must not run for ghost points.
1097    CALL module_interface_exchange_horiz()
1098!
1099!-- u-velocity component
1100    CALL cpu_log( log_point(5), 'u-equation', 'start' )
1101
1102    !$ACC KERNELS PRESENT(tend)
1103    tend = 0.0_wp
1104    !$ACC END KERNELS
1105    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1106       IF ( ws_scheme_mom )  THEN
1107          CALL advec_u_ws
1108       ELSE
1109          CALL advec_u_pw
1110       ENDIF
1111    ELSE
1112       CALL advec_u_up
1113    ENDIF
1114    CALL diffusion_u
1115    CALL coriolis( 1 )
1116    IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
1117       CALL buoyancy( pt, 1 )
1118    ENDIF
1119
1120!
1121!-- Drag by plant canopy
1122    IF ( plant_canopy )  CALL pcm_tendency( 1 )
1123
1124!
1125!-- External pressure gradient
1126    IF ( dp_external )  THEN
1127       DO  i = nxlu, nxr
1128          DO  j = nys, nyn
1129             DO  k = dp_level_ind_b+1, nzt
1130                tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
1131             ENDDO
1132          ENDDO
1133       ENDDO
1134    ENDIF
1135
1136!
1137!-- Nudging
1138    IF ( nudging )  CALL nudge( simulated_time, 'u' )
1139
1140!
1141!-- Effect of Stokes drift (in ocean mode only)
1142    IF ( stokes_force )  CALL stokes_drift_terms( 1 )
1143
1144    CALL module_interface_actions( 'u-tendency' )
1145
1146!
1147!-- Prognostic equation for u-velocity component
1148    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1149    !$ACC PRESENT(u, tend, tu_m, u_init, rdf, wall_flags_0) &
1150    !$ACC PRESENT(tsc(2:5)) &
1151    !$ACC PRESENT(u_p)
1152    DO  i = nxlu, nxr
1153       DO  j = nys, nyn
1154          DO  k = nzb+1, nzt
1155             u_p(k,j,i) = u(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +          &
1156                                                 tsc(3) * tu_m(k,j,i) )          &
1157                                               - tsc(5) * rdf(k) *               &
1158                                                        ( u(k,j,i) - u_init(k) ) &
1159                                     ) * MERGE( 1.0_wp, 0.0_wp,                  &
1160                                                 BTEST( wall_flags_0(k,j,i), 1 ) &
1161                                              )
1162          ENDDO
1163       ENDDO
1164    ENDDO
1165
1166!
1167!-- Add turbulence generated by wave breaking (in ocean mode only)
1168    IF ( wave_breaking  .AND.                                                  &
1169         intermediate_timestep_count == intermediate_timestep_count_max )      &
1170    THEN
1171       CALL wave_breaking_term( 1 )
1172    ENDIF
1173
1174!
1175!-- Calculate tendencies for the next Runge-Kutta step
1176    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1177       IF ( intermediate_timestep_count == 1 )  THEN
1178          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1179          !$ACC PRESENT(tend, tu_m)
1180          DO  i = nxlu, nxr
1181             DO  j = nys, nyn
1182                DO  k = nzb+1, nzt
1183                   tu_m(k,j,i) = tend(k,j,i)
1184                ENDDO
1185             ENDDO
1186          ENDDO
1187       ELSEIF ( intermediate_timestep_count < &
1188                intermediate_timestep_count_max )  THEN
1189          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1190          !$ACC PRESENT(tend, tu_m)
1191          DO  i = nxlu, nxr
1192             DO  j = nys, nyn
1193                DO  k = nzb+1, nzt
1194                   tu_m(k,j,i) =    -9.5625_wp * tend(k,j,i)                   &
1195                                   + 5.3125_wp * tu_m(k,j,i)
1196                ENDDO
1197             ENDDO
1198          ENDDO
1199       ENDIF
1200    ENDIF
1201
1202    CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1203
1204!
1205!-- v-velocity component
1206    CALL cpu_log( log_point(6), 'v-equation', 'start' )
1207
1208    !$ACC KERNELS PRESENT(tend)
1209    tend = 0.0_wp
1210    !$ACC END KERNELS
1211    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1212       IF ( ws_scheme_mom )  THEN
1213          CALL advec_v_ws
1214       ELSE
1215          CALL advec_v_pw
1216       END IF
1217    ELSE
1218       CALL advec_v_up
1219    ENDIF
1220    CALL diffusion_v
1221    CALL coriolis( 2 )
1222
1223!
1224!-- Drag by plant canopy
1225    IF ( plant_canopy )  CALL pcm_tendency( 2 )
1226
1227!
1228!-- External pressure gradient
1229    IF ( dp_external )  THEN
1230       DO  i = nxl, nxr
1231          DO  j = nysv, nyn
1232             DO  k = dp_level_ind_b+1, nzt
1233                tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1234             ENDDO
1235          ENDDO
1236       ENDDO
1237    ENDIF
1238
1239!
1240!-- Nudging
1241    IF ( nudging )  CALL nudge( simulated_time, 'v' )
1242
1243!
1244!-- Effect of Stokes drift (in ocean mode only)
1245    IF ( stokes_force )  CALL stokes_drift_terms( 2 )
1246
1247    CALL module_interface_actions( 'v-tendency' )
1248
1249!
1250!-- Prognostic equation for v-velocity component
1251    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1252    !$ACC PRESENT(v, tend, tv_m, v_init, rdf, wall_flags_0) &
1253    !$ACC PRESENT(tsc(2:5)) &
1254    !$ACC PRESENT(v_p)
1255    DO  i = nxl, nxr
1256       DO  j = nysv, nyn
1257          DO  k = nzb+1, nzt
1258             v_p(k,j,i) = v(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1259                                                 tsc(3) * tv_m(k,j,i) )        &
1260                                               - tsc(5) * rdf(k) *             &
1261                                                      ( v(k,j,i) - v_init(k) ) &
1262                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1263                                                BTEST( wall_flags_0(k,j,i), 2 )&
1264                                              )
1265          ENDDO
1266       ENDDO
1267    ENDDO
1268
1269!
1270!-- Add turbulence generated by wave breaking (in ocean mode only)
1271    IF ( wave_breaking  .AND.                                                  &
1272         intermediate_timestep_count == intermediate_timestep_count_max )      &
1273    THEN
1274       CALL wave_breaking_term( 2 )
1275    ENDIF
1276
1277!
1278!-- Calculate tendencies for the next Runge-Kutta step
1279    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1280       IF ( intermediate_timestep_count == 1 )  THEN
1281          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1282          !$ACC PRESENT(tend, tv_m)
1283          DO  i = nxl, nxr
1284             DO  j = nysv, nyn
1285                DO  k = nzb+1, nzt
1286                   tv_m(k,j,i) = tend(k,j,i)
1287                ENDDO
1288             ENDDO
1289          ENDDO
1290       ELSEIF ( intermediate_timestep_count < &
1291                intermediate_timestep_count_max )  THEN
1292          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1293          !$ACC PRESENT(tend, tv_m)
1294          DO  i = nxl, nxr
1295             DO  j = nysv, nyn
1296                DO  k = nzb+1, nzt
1297                   tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1298                                  + 5.3125_wp * tv_m(k,j,i)
1299                ENDDO
1300             ENDDO
1301          ENDDO
1302       ENDIF
1303    ENDIF
1304
1305    CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1306
1307!
1308!-- w-velocity component
1309    CALL cpu_log( log_point(7), 'w-equation', 'start' )
1310
1311    !$ACC KERNELS PRESENT(tend)
1312    tend = 0.0_wp
1313    !$ACC END KERNELS
1314    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1315       IF ( ws_scheme_mom )  THEN
1316          CALL advec_w_ws
1317       ELSE
1318          CALL advec_w_pw
1319       ENDIF
1320    ELSE
1321       CALL advec_w_up
1322    ENDIF
1323    CALL diffusion_w
1324    CALL coriolis( 3 )
1325
1326    IF ( .NOT. neutral )  THEN
1327       IF ( ocean_mode )  THEN
1328          CALL buoyancy( rho_ocean, 3 )
1329       ELSE
1330          IF ( .NOT. humidity )  THEN
1331             CALL buoyancy( pt, 3 )
1332          ELSE
1333             CALL buoyancy( vpt, 3 )
1334          ENDIF
1335       ENDIF
1336    ENDIF
1337
1338!
1339!-- Drag by plant canopy
1340    IF ( plant_canopy )  CALL pcm_tendency( 3 )
1341
1342!
1343!-- Effect of Stokes drift (in ocean mode only)
1344    IF ( stokes_force )  CALL stokes_drift_terms( 3 )
1345
1346    CALL module_interface_actions( 'w-tendency' )
1347
1348!
1349!-- Prognostic equation for w-velocity component
1350    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1351    !$ACC PRESENT(w, tend, tw_m, v_init, rdf, wall_flags_0) &
1352    !$ACC PRESENT(tsc(2:5)) &
1353    !$ACC PRESENT(w_p)
1354    DO  i = nxl, nxr
1355       DO  j = nys, nyn
1356          DO  k = nzb+1, nzt-1
1357             w_p(k,j,i) = w(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1358                                                 tsc(3) * tw_m(k,j,i) )        &
1359                                               - tsc(5) * rdf(k) * w(k,j,i)    &
1360                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1361                                                BTEST( wall_flags_0(k,j,i), 3 )&
1362                                              )
1363          ENDDO
1364       ENDDO
1365    ENDDO
1366
1367!
1368!-- Calculate tendencies for the next Runge-Kutta step
1369    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1370       IF ( intermediate_timestep_count == 1 )  THEN
1371          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1372          !$ACC PRESENT(tend, tw_m)
1373          DO  i = nxl, nxr
1374             DO  j = nys, nyn
1375                DO  k = nzb+1, nzt-1
1376                   tw_m(k,j,i) = tend(k,j,i)
1377                ENDDO
1378             ENDDO
1379          ENDDO
1380       ELSEIF ( intermediate_timestep_count < &
1381                intermediate_timestep_count_max )  THEN
1382          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1383          !$ACC PRESENT(tend, tw_m)
1384          DO  i = nxl, nxr
1385             DO  j = nys, nyn
1386                DO  k = nzb+1, nzt-1
1387                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1388                                  + 5.3125_wp * tw_m(k,j,i)
1389                ENDDO
1390             ENDDO
1391          ENDDO
1392       ENDIF
1393    ENDIF
1394
1395    CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1396
1397
1398!
1399!-- If required, compute prognostic equation for potential temperature
1400    IF ( .NOT. neutral )  THEN
1401
1402       CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1403
1404!
1405!--    pt-tendency terms with communication
1406       sbt = tsc(2)
1407       IF ( scalar_advec == 'bc-scheme' )  THEN
1408
1409          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1410!
1411!--          Bott-Chlond scheme always uses Euler time step. Thus:
1412             sbt = 1.0_wp
1413          ENDIF
1414          tend = 0.0_wp
1415          CALL advec_s_bc( pt, 'pt' )
1416
1417       ENDIF
1418
1419!
1420!--    pt-tendency terms with no communication
1421       IF ( scalar_advec /= 'bc-scheme' )  THEN
1422          !$ACC KERNELS PRESENT(tend)
1423          tend = 0.0_wp
1424          !$ACC END KERNELS
1425          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1426             IF ( ws_scheme_sca )  THEN
1427                CALL advec_s_ws( pt, 'pt' )
1428             ELSE
1429                CALL advec_s_pw( pt )
1430             ENDIF
1431          ELSE
1432             CALL advec_s_up( pt )
1433          ENDIF
1434       ENDIF
1435
1436       CALL diffusion_s( pt,                                                   &
1437                         surf_def_h(0)%shf, surf_def_h(1)%shf,                 &
1438                         surf_def_h(2)%shf,                                    &
1439                         surf_lsm_h%shf,    surf_usm_h%shf,                    &
1440                         surf_def_v(0)%shf, surf_def_v(1)%shf,                 &
1441                         surf_def_v(2)%shf, surf_def_v(3)%shf,                 &
1442                         surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,                 &
1443                         surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,                 &
1444                         surf_usm_v(0)%shf, surf_usm_v(1)%shf,                 &
1445                         surf_usm_v(2)%shf, surf_usm_v(3)%shf )
1446
1447!
1448!--    Consideration of heat sources within the plant canopy
1449       IF ( plant_canopy  .AND.                                          &
1450            (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
1451          CALL pcm_tendency( 4 )
1452       ENDIF
1453
1454!
1455!--    Large scale advection
1456       IF ( large_scale_forcing )  THEN
1457          CALL ls_advec( simulated_time, 'pt' )
1458       ENDIF
1459
1460!
1461!--    Nudging
1462       IF ( nudging )  CALL nudge( simulated_time, 'pt' )
1463
1464!
1465!--    If required compute influence of large-scale subsidence/ascent
1466       IF ( large_scale_subsidence  .AND.                                      &
1467            .NOT. use_subsidence_tendencies )  THEN
1468          CALL subsidence( tend, pt, pt_init, 2 )
1469       ENDIF
1470
1471#if defined( __rrtmg )
1472!
1473!--    If required, add tendency due to radiative heating/cooling
1474       IF ( radiation  .AND.                                                   &
1475            simulated_time > skip_time_do_radiation )  THEN
1476            CALL radiation_tendency ( tend )
1477       ENDIF
1478#endif
1479
1480       CALL module_interface_actions( 'pt-tendency' )
1481
1482!
1483!--    Prognostic equation for potential temperature
1484       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1485       !$ACC PRESENT(pt, tend, tpt_m, wall_flags_0) &
1486       !$ACC PRESENT(pt_init, rdf_sc, ptdf_x, ptdf_y) &
1487       !$ACC PRESENT(tsc(3:5)) &
1488       !$ACC PRESENT(pt_p)
1489       DO  i = nxl, nxr
1490          DO  j = nys, nyn
1491             DO  k = nzb+1, nzt
1492                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +      &
1493                                                   tsc(3) * tpt_m(k,j,i) )     &
1494                                                 - tsc(5) *                    &
1495                                                   ( pt(k,j,i) - pt_init(k) ) *&
1496                                          ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )&
1497                                          )                                    &
1498                                   * MERGE( 1.0_wp, 0.0_wp,                    &
1499                                             BTEST( wall_flags_0(k,j,i), 0 )   &
1500                                          )
1501             ENDDO
1502          ENDDO
1503       ENDDO
1504!
1505!--    Calculate tendencies for the next Runge-Kutta step
1506       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1507          IF ( intermediate_timestep_count == 1 )  THEN
1508             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1509             !$ACC PRESENT(tend, tpt_m)
1510             DO  i = nxl, nxr
1511                DO  j = nys, nyn
1512                   DO  k = nzb+1, nzt
1513                      tpt_m(k,j,i) = tend(k,j,i)
1514                   ENDDO
1515                ENDDO
1516             ENDDO
1517          ELSEIF ( intermediate_timestep_count < &
1518                   intermediate_timestep_count_max )  THEN
1519             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1520             !$ACC PRESENT(tend, tpt_m)
1521             DO  i = nxl, nxr
1522                DO  j = nys, nyn
1523                   DO  k = nzb+1, nzt
1524                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
1525                                        5.3125_wp * tpt_m(k,j,i)
1526                   ENDDO
1527                ENDDO
1528             ENDDO
1529          ENDIF
1530       ENDIF
1531
1532       CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1533
1534    ENDIF
1535
1536!
1537!-- If required, compute prognostic equation for total water content
1538    IF ( humidity )  THEN
1539
1540       CALL cpu_log( log_point(29), 'q-equation', 'start' )
1541
1542!
1543!--    Scalar/q-tendency terms with communication
1544       sbt = tsc(2)
1545       IF ( scalar_advec == 'bc-scheme' )  THEN
1546
1547          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1548!
1549!--          Bott-Chlond scheme always uses Euler time step. Thus:
1550             sbt = 1.0_wp
1551          ENDIF
1552          tend = 0.0_wp
1553          CALL advec_s_bc( q, 'q' )
1554
1555       ENDIF
1556
1557!
1558!--    Scalar/q-tendency terms with no communication
1559       IF ( scalar_advec /= 'bc-scheme' )  THEN
1560          tend = 0.0_wp
1561          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1562             IF ( ws_scheme_sca )  THEN
1563                CALL advec_s_ws( q, 'q' )
1564             ELSE
1565                CALL advec_s_pw( q )
1566             ENDIF
1567          ELSE
1568             CALL advec_s_up( q )
1569          ENDIF
1570       ENDIF
1571
1572       CALL diffusion_s( q,                                                    &
1573                         surf_def_h(0)%qsws, surf_def_h(1)%qsws,               &
1574                         surf_def_h(2)%qsws,                                   &
1575                         surf_lsm_h%qsws,    surf_usm_h%qsws,                  &
1576                         surf_def_v(0)%qsws, surf_def_v(1)%qsws,               &
1577                         surf_def_v(2)%qsws, surf_def_v(3)%qsws,               &
1578                         surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,               &
1579                         surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,               &
1580                         surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,               &
1581                         surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
1582
1583!
1584!--    Sink or source of humidity due to canopy elements
1585       IF ( plant_canopy ) CALL pcm_tendency( 5 )
1586
1587!
1588!--    Large scale advection
1589       IF ( large_scale_forcing )  THEN
1590          CALL ls_advec( simulated_time, 'q' )
1591       ENDIF
1592
1593!
1594!--    Nudging
1595       IF ( nudging )  CALL nudge( simulated_time, 'q' )
1596
1597!
1598!--    If required compute influence of large-scale subsidence/ascent
1599       IF ( large_scale_subsidence  .AND.                                      &
1600            .NOT. use_subsidence_tendencies )  THEN
1601         CALL subsidence( tend, q, q_init, 3 )
1602       ENDIF
1603
1604       CALL module_interface_actions( 'q-tendency' )
1605
1606!
1607!--    Prognostic equation for total water content
1608       DO  i = nxl, nxr
1609          DO  j = nys, nyn
1610             DO  k = nzb+1, nzt
1611                q_p(k,j,i) = q(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
1612                                                 tsc(3) * tq_m(k,j,i) )        &
1613                                               - tsc(5) * rdf_sc(k) *          &
1614                                                      ( q(k,j,i) - q_init(k) ) &
1615                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
1616                                                   BTEST( wall_flags_0(k,j,i), 0 ) &
1617                                                 )
1618                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
1619             ENDDO
1620          ENDDO
1621       ENDDO
1622
1623!
1624!--    Calculate tendencies for the next Runge-Kutta step
1625       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1626          IF ( intermediate_timestep_count == 1 )  THEN
1627             DO  i = nxl, nxr
1628                DO  j = nys, nyn
1629                   DO  k = nzb+1, nzt
1630                      tq_m(k,j,i) = tend(k,j,i)
1631                   ENDDO
1632                ENDDO
1633             ENDDO
1634          ELSEIF ( intermediate_timestep_count < &
1635                   intermediate_timestep_count_max )  THEN
1636             DO  i = nxl, nxr
1637                DO  j = nys, nyn
1638                   DO  k = nzb+1, nzt
1639                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
1640                                     + 5.3125_wp * tq_m(k,j,i)
1641                   ENDDO
1642                ENDDO
1643             ENDDO
1644          ENDIF
1645       ENDIF
1646
1647       CALL cpu_log( log_point(29), 'q-equation', 'stop' )
1648
1649    ENDIF
1650!
1651!-- If required, compute prognostic equation for scalar
1652    IF ( passive_scalar )  THEN
1653
1654       CALL cpu_log( log_point(66), 's-equation', 'start' )
1655
1656!
1657!--    Scalar/q-tendency terms with communication
1658       sbt = tsc(2)
1659       IF ( scalar_advec == 'bc-scheme' )  THEN
1660
1661          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1662!
1663!--          Bott-Chlond scheme always uses Euler time step. Thus:
1664             sbt = 1.0_wp
1665          ENDIF
1666          tend = 0.0_wp
1667          CALL advec_s_bc( s, 's' )
1668
1669       ENDIF
1670
1671!
1672!--    Scalar/q-tendency terms with no communication
1673       IF ( scalar_advec /= 'bc-scheme' )  THEN
1674          tend = 0.0_wp
1675          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1676             IF ( ws_scheme_sca )  THEN
1677                CALL advec_s_ws( s, 's' )
1678             ELSE
1679                CALL advec_s_pw( s )
1680             ENDIF
1681          ELSE
1682             CALL advec_s_up( s )
1683          ENDIF
1684       ENDIF
1685
1686       CALL diffusion_s( s,                                                    &
1687                         surf_def_h(0)%ssws, surf_def_h(1)%ssws,               &
1688                         surf_def_h(2)%ssws,                                   &
1689                         surf_lsm_h%ssws,    surf_usm_h%ssws,                  &
1690                         surf_def_v(0)%ssws, surf_def_v(1)%ssws,               &
1691                         surf_def_v(2)%ssws, surf_def_v(3)%ssws,               &
1692                         surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,               &
1693                         surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,               &
1694                         surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,               &
1695                         surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1696
1697!
1698!--    Sink or source of humidity due to canopy elements
1699       IF ( plant_canopy ) CALL pcm_tendency( 7 )
1700
1701!
1702!--    Large scale advection. Not implemented for scalars so far.
1703!        IF ( large_scale_forcing )  THEN
1704!           CALL ls_advec( simulated_time, 'q' )
1705!        ENDIF
1706
1707!
1708!--    Nudging. Not implemented for scalars so far.
1709!        IF ( nudging )  CALL nudge( simulated_time, 'q' )
1710
1711!
1712!--    If required compute influence of large-scale subsidence/ascent.
1713!--    Not implemented for scalars so far.
1714       IF ( large_scale_subsidence  .AND.                                      &
1715            .NOT. use_subsidence_tendencies  .AND.                             &
1716            .NOT. large_scale_forcing )  THEN
1717         CALL subsidence( tend, s, s_init, 3 )
1718       ENDIF
1719
1720       CALL module_interface_actions( 's-tendency' )
1721
1722!
1723!--    Prognostic equation for total water content
1724       DO  i = nxl, nxr
1725          DO  j = nys, nyn
1726             DO  k = nzb+1, nzt
1727                s_p(k,j,i) = s(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
1728                                                 tsc(3) * ts_m(k,j,i) )        &
1729                                               - tsc(5) * rdf_sc(k) *          &
1730                                                    ( s(k,j,i) - s_init(k) )   &
1731                                        )                                      &
1732                                   * MERGE( 1.0_wp, 0.0_wp,                    &
1733                                             BTEST( wall_flags_0(k,j,i), 0 )   &
1734                                          )
1735                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
1736             ENDDO
1737          ENDDO
1738       ENDDO
1739
1740!
1741!--    Calculate tendencies for the next Runge-Kutta step
1742       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1743          IF ( intermediate_timestep_count == 1 )  THEN
1744             DO  i = nxl, nxr
1745                DO  j = nys, nyn
1746                   DO  k = nzb+1, nzt
1747                      ts_m(k,j,i) = tend(k,j,i)
1748                   ENDDO
1749                ENDDO
1750             ENDDO
1751          ELSEIF ( intermediate_timestep_count < &
1752                   intermediate_timestep_count_max )  THEN
1753             DO  i = nxl, nxr
1754                DO  j = nys, nyn
1755                   DO  k = nzb+1, nzt
1756                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
1757                                     + 5.3125_wp * ts_m(k,j,i)
1758                   ENDDO
1759                ENDDO
1760             ENDDO
1761          ENDIF
1762       ENDIF
1763
1764       CALL cpu_log( log_point(66), 's-equation', 'stop' )
1765
1766    ENDIF
1767!
1768!-- Calculate prognostic equations for all other modules
1769    CALL module_interface_prognostic_equations()
1770
1771    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_vector', 'end' )
1772
1773 END SUBROUTINE prognostic_equations_vector
1774
1775
1776 END MODULE prognostic_equations_mod
Note: See TracBrowser for help on using the repository browser.