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

Last change on this file since 4028 was 3987, checked in by kanani, 5 years ago

clean up location, debug and error messages

  • Property svn:keywords set to Id
File size: 64.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 3987 2019-05-22 09:52:13Z schwenkel $
27! Introduce alternative switch for debug output during timestepping
28!
29! 3956 2019-05-07 12:32:52Z monakurppa
30! Removed salsa calls.
31!
32! 3931 2019-04-24 16:34:28Z schwenkel
33! Correct/complete module_interface introduction for chemistry model
34!
35! 3899 2019-04-16 14:05:27Z monakurppa
36! Corrections in the OpenMP version of salsa
37!
38! 3887 2019 -04-12 08:47:41Z schwenkel
39! Implicit Bugfix for chemistry model, loop for non_transport_physics over
40! ghost points is avoided. Instead introducing module_interface_exchange_horiz.
41!
42! 3885 2019-04-11 11:29:34Z kanani
43! Changes related to global restructuring of location messages and introduction
44! of additional debug messages
45!
46! 3881 2019-04-10 09:31:22Z suehring
47! Bugfix in OpenMP directive
48!
49! 3880 2019-04-08 21:43:02Z knoop
50! Moved wtm_tendencies to module_interface_actions
51!
52! 3874 2019-04-08 16:53:48Z knoop
53! Added non_transport_physics module interfaces and moved bcm code into it
54!
55! 3872 2019-04-08 15:03:06Z knoop
56! Moving prognostic equations of bcm into bulk_cloud_model_mod
57!
58! 3864 2019-04-05 09:01:56Z monakurppa
59! Modifications made for salsa:
60! - salsa_prognostic_equations moved to salsa_mod (and the call to
61!   module_interface_mod)
62! - Renamed nbins --> nbins_aerosol, ncc_tot --> ncomponents_mass and
63!   ngast --> ngases_salsa and loop indices b, c and sg to ib, ic and ig
64!
65! 3840 2019-03-29 10:35:52Z knoop
66! added USE chem_gasphase_mod for nspec, nspec and spc_names
67!
68! 3820 2019-03-27 11:53:41Z forkel
69! renamed do_depo to deposition_dry (ecc)
70!
71! 3797 2019-03-15 11:15:38Z forkel
72! Call chem_integegrate in OpenMP loop   (ketelsen)
73!
74!
75! 3771 2019-02-28 12:19:33Z raasch
76! preprocessor directivs fro rrtmg added
77!
78! 3761 2019-02-25 15:31:42Z raasch
79! unused variable removed
80!
81! 3719 2019-02-06 13:10:18Z kanani
82! Cleaned up chemistry cpu measurements
83!
84! 3684 2019-01-20 20:20:58Z knoop
85! OpenACC port for SPEC
86!
87! 3589 2018-11-30 15:09:51Z suehring
88! Move the control parameter "salsa" from salsa_mod to control_parameters
89! (M. Kurppa)
90!
91! 3582 2018-11-29 19:16:36Z suehring
92! Implementation of a new aerosol module salsa.
93! Remove cpu-logs from i,j loop in cache version.
94!
95! 3458 2018-10-30 14:51:23Z kanani
96! remove duplicate USE chem_modules
97! from chemistry branch r3443, banzhafs, basit:
98! chem_depo call introduced
99! code added for decycling chemistry
100!
101! 3386 2018-10-19 16:28:22Z gronemeier
102! Renamed tcm_prognostic to tcm_prognostic_equations
103!
104! 3355 2018-10-16 14:03:34Z knoop
105! (from branch resler)
106! Fix for chemistry call
107!
108! 3302 2018-10-03 02:39:40Z raasch
109! Stokes drift + wave breaking term added
110!
111! 3298 2018-10-02 12:21:11Z kanani
112! Code added for decycling chemistry (basit)
113!
114! 3294 2018-10-01 02:37:10Z raasch
115! changes concerning modularization of ocean option
116!
117! 3274 2018-09-24 15:42:55Z knoop
118! Modularization of all bulk cloud physics code components
119!
120! 3241 2018-09-12 15:02:00Z raasch
121! omp_get_thread_num now declared in omp directive
122!
123! 3183 2018-07-27 14:25:55Z suehring
124! Remove unused variables from USE statements
125!
126! 3182 2018-07-27 13:36:03Z suehring
127! Revise recent bugfix for nesting
128!
129! 3021 2018-05-16 08:14:20Z maronga
130! Bugfix in IF clause for nesting
131!
132! 3014 2018-05-09 08:42:38Z maronga
133! Fixed a bug in the IF condition to call pcm_tendency in case of
134! potential temperature
135!
136! 2815 2018-02-19 11:29:57Z kanani
137! Rename chem_tendency to chem_prognostic_equations,
138! implement vector version for air chemistry
139!
140! 2766 2018-01-22 17:17:47Z kanani
141! Removed preprocessor directive __chem
142!
143! 2746 2018-01-15 12:06:04Z suehring
144! Move flag plant canopy to modules
145!
146! 2719 2018-01-02 09:02:06Z maronga
147! Bugfix for last change.
148!
149! 2718 2018-01-02 08:49:38Z maronga
150! Corrected "Former revisions" section
151!
152! 2696 2017-12-14 17:12:51Z kanani
153! - Change in file header (GPL part)
154! - Moved TKE equation to tcm_prognostic (TG)
155! - Added switch for chemical reactions (RF, FK)
156! - Implementation of chemistry module (RF, BK, FK)
157!
158! 2563 2017-10-19 15:36:10Z Giersch
159! Variable wind_turbine moved to module control_parameters
160!
161! 2320 2017-07-21 12:47:43Z suehring
162! Modularize large-scale forcing and nudging
163!
164! 2292 2017-06-20 09:51:42Z schwenkel
165! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
166! includes two more prognostic equations for cloud drop concentration (nc) 
167! and cloud water content (qc).
168!
169! 2261 2017-06-08 14:25:57Z raasch
170! bugfix for r2232: openmp directives removed
171!
172! 2233 2017-05-30 18:08:54Z suehring
173!
174! 2232 2017-05-30 17:47:52Z suehring
175! Adjutst to new surface-type structure. Remove call for usm_wall_heat_flux,
176! which is realized directly in diffusion_s now.
177!
178! 2192 2017-03-22 04:14:10Z raasch
179! Bugfix for misplaced and missing openMP directives from r2155
180!
181! 2155 2017-02-21 09:57:40Z hoffmann
182! Bugfix in the calculation of microphysical quantities on ghost points.
183!
184! 2118 2017-01-17 16:38:49Z raasch
185! OpenACC version of subroutine removed
186!
187! 2031 2016-10-21 15:11:58Z knoop
188! renamed variable rho to rho_ocean
189!
190! 2011 2016-09-19 17:29:57Z kanani
191! Flag urban_surface is now defined in module control_parameters.
192!
193! 2007 2016-08-24 15:47:17Z kanani
194! Added pt tendency calculation based on energy balance at urban surfaces
195! (new urban surface model)
196!
197! 2000 2016-08-20 18:09:15Z knoop
198! Forced header and separation lines into 80 columns
199!
200! 1976 2016-07-27 13:28:04Z maronga
201! Simplied calls to radiation model
202!
203! 1960 2016-07-12 16:34:24Z suehring
204! Separate humidity and passive scalar
205!
206! 1914 2016-05-26 14:44:07Z witha
207! Added calls for wind turbine model
208!
209! 1873 2016-04-18 14:50:06Z maronga
210! Module renamed (removed _mod)
211!
212! 1850 2016-04-08 13:29:27Z maronga
213! Module renamed
214!
215! 1826 2016-04-07 12:01:39Z maronga
216! Renamed canopy model calls.
217!
218! 1822 2016-04-07 07:49:42Z hoffmann
219! Kessler microphysics scheme moved to microphysics.
220!
221! 1757 2016-02-22 15:49:32Z maronga
222!
223! 1691 2015-10-26 16:17:44Z maronga
224! Added optional model spin-up without radiation / land surface model calls.
225! Formatting corrections.
226!
227! 1682 2015-10-07 23:56:08Z knoop
228! Code annotations made doxygen readable
229!
230! 1585 2015-04-30 07:05:52Z maronga
231! Added call for temperature tendency calculation due to radiative flux divergence
232!
233! 1517 2015-01-07 19:12:25Z hoffmann
234! advec_s_bc_mod addded, since advec_s_bc is now a module
235!
236! 1484 2014-10-21 10:53:05Z kanani
237! Changes due to new module structure of the plant canopy model:
238! parameters cthf and plant_canopy moved to module plant_canopy_model_mod.
239! Removed double-listing of use_upstream_for_tke in ONLY-list of module
240! control_parameters
241!
242! 1409 2014-05-23 12:11:32Z suehring
243! Bugfix: i_omp_start changed for advec_u_ws at left inflow and outflow boundary.
244! This ensures that left-hand side fluxes are also calculated for nxl in that
245! case, even though the solution at nxl is overwritten in boundary_conds()
246!
247! 1398 2014-05-07 11:15:00Z heinze
248! Rayleigh-damping for horizontal velocity components changed: instead of damping
249! against ug and vg, damping against u_init and v_init is used to allow for a
250! homogenized treatment in case of nudging
251!
252! 1380 2014-04-28 12:40:45Z heinze
253! Change order of calls for scalar prognostic quantities:
254! ls_advec -> nudging -> subsidence since initial profiles
255!
256! 1374 2014-04-25 12:55:07Z raasch
257! missing variables added to ONLY lists
258!
259! 1365 2014-04-22 15:03:56Z boeske
260! Calls of ls_advec for large scale advection added,
261! subroutine subsidence is only called if use_subsidence_tendencies = .F.,
262! new argument ls_index added to the calls of subsidence
263! +ls_index
264!
265! 1361 2014-04-16 15:17:48Z hoffmann
266! Two-moment microphysics moved to the start of prognostic equations. This makes
267! the 3d arrays for tend_q, tend_qr, tend_pt and tend_pt redundant.
268! Additionally, it is allowed to call the microphysics just once during the time
269! step (not at each sub-time step).
270!
271! Two-moment cloud physics added for vector and accelerator optimization.
272!
273! 1353 2014-04-08 15:21:23Z heinze
274! REAL constants provided with KIND-attribute
275!
276! 1337 2014-03-25 15:11:48Z heinze
277! Bugfix: REAL constants provided with KIND-attribute
278!
279! 1332 2014-03-25 11:59:43Z suehring
280! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke
281!
282! 1330 2014-03-24 17:29:32Z suehring
283! In case of SGS-particle velocity advection of TKE is also allowed with
284! dissipative 5th-order scheme.
285!
286! 1320 2014-03-20 08:40:49Z raasch
287! ONLY-attribute added to USE-statements,
288! kind-parameters added to all INTEGER and REAL declaration statements,
289! kinds are defined in new module kinds,
290! old module precision_kind is removed,
291! revision history before 2012 removed,
292! comment fields (!:) to be used for variable explanations added to
293! all variable declaration statements
294!
295! 1318 2014-03-17 13:35:16Z raasch
296! module interfaces removed
297!
298! 1257 2013-11-08 15:18:40Z raasch
299! openacc loop vector clauses removed, independent clauses added
300!
301! 1246 2013-11-01 08:59:45Z heinze
302! enable nudging also for accelerator version
303!
304! 1241 2013-10-30 11:36:58Z heinze
305! usage of nudging enabled (so far not implemented for accelerator version)
306!
307! 1179 2013-06-14 05:57:58Z raasch
308! two arguments removed from routine buoyancy, ref_state updated on device
309!
310! 1128 2013-04-12 06:19:32Z raasch
311! those parts requiring global communication moved to time_integration,
312! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
313! j_north
314!
315! 1115 2013-03-26 18:16:16Z hoffmann
316! optimized cloud physics: calculation of microphysical tendencies transfered
317! to microphysics.f90; qr and nr are only calculated if precipitation is required
318!
319! 1111 2013-03-08 23:54:10Z raasch
320! update directives for prognostic quantities removed
321!
322! 1106 2013-03-04 05:31:38Z raasch
323! small changes in code formatting
324!
325! 1092 2013-02-02 11:24:22Z raasch
326! unused variables removed
327!
328! 1053 2012-11-13 17:11:03Z hoffmann
329! implementation of two new prognostic equations for rain drop concentration (nr)
330! and rain water content (qr)
331!
332! currently, only available for cache loop optimization
333!
334! 1036 2012-10-22 13:43:42Z raasch
335! code put under GPL (PALM 3.9)
336!
337! 1019 2012-09-28 06:46:45Z raasch
338! non-optimized version of prognostic_equations removed
339!
340! 1015 2012-09-27 09:23:24Z raasch
341! new branch prognostic_equations_acc
342! OpenACC statements added + code changes required for GPU optimization
343!
344! 1001 2012-09-13 14:08:46Z raasch
345! all actions concerning leapfrog- and upstream-spline-scheme removed
346!
347! 978 2012-08-09 08:28:32Z fricke
348! km_damp_x and km_damp_y removed in calls of diffusion_u and diffusion_v
349! add ptdf_x, ptdf_y for damping the potential temperature at the inflow
350! boundary in case of non-cyclic lateral boundaries
351! Bugfix: first thread index changes for WS-scheme at the inflow
352!
353! 940 2012-07-09 14:31:00Z raasch
354! temperature equation can be switched off
355!
356! Revision 1.1  2000/04/13 14:56:27  schroeter
357! Initial revision
358!
359!
360! Description:
361! ------------
362!> Solving the prognostic equations.
363!------------------------------------------------------------------------------!
364 MODULE prognostic_equations_mod
365
366    USE advec_s_bc_mod,                                                        &
367        ONLY:  advec_s_bc
368
369    USE advec_s_pw_mod,                                                        &
370        ONLY:  advec_s_pw
371
372    USE advec_s_up_mod,                                                        &
373        ONLY:  advec_s_up
374
375    USE advec_u_pw_mod,                                                        &
376        ONLY:  advec_u_pw
377
378    USE advec_u_up_mod,                                                        &
379        ONLY:  advec_u_up
380
381    USE advec_v_pw_mod,                                                        &
382        ONLY:  advec_v_pw
383
384    USE advec_v_up_mod,                                                        &
385        ONLY:  advec_v_up
386
387    USE advec_w_pw_mod,                                                        &
388        ONLY:  advec_w_pw
389
390    USE advec_w_up_mod,                                                        &
391        ONLY:  advec_w_up
392
393    USE advec_ws,                                                              &
394        ONLY:  advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws
395
396    USE arrays_3d,                                                             &
397        ONLY:  diss_l_e, diss_l_pt, diss_l_q,                                  &
398               diss_l_s, diss_l_sa, diss_s_e,                                  &
399               diss_s_pt, diss_s_q, diss_s_s, diss_s_sa,                       &
400               e, e_p, flux_s_e, flux_s_pt, flux_s_q,                          &
401               flux_s_s, flux_s_sa, flux_l_e,                                  &
402               flux_l_pt, flux_l_q, flux_l_s,                                  &
403               flux_l_sa, pt, ptdf_x, ptdf_y, pt_init,                         &
404               pt_p, prho, q, q_init, q_p, rdf, rdf_sc,                        &
405               ref_state, rho_ocean, s, s_init, s_p, tend, te_m,               &
406               tpt_m, tq_m, ts_m, tu_m, tv_m, tw_m, u,                         &
407               ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
408
409    USE buoyancy_mod,                                                          &
410        ONLY:  buoyancy
411
412    USE control_parameters,                                                    &
413        ONLY:  constant_diffusion,                                             &
414               debug_output_timestep,                                          &
415               dp_external, dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d,    &
416               humidity, intermediate_timestep_count,                          &
417               intermediate_timestep_count_max, large_scale_forcing,           &
418               large_scale_subsidence, neutral, nudging,                       &
419               ocean_mode, passive_scalar, plant_canopy, pt_reference,         &
420               scalar_advec, scalar_advec, simulated_time, sloping_surface,    &
421               timestep_scheme, tsc, use_subsidence_tendencies,                &
422               use_upstream_for_tke, wind_turbine, ws_scheme_mom,              & 
423               ws_scheme_sca, urban_surface, land_surface,                     &
424               time_since_reference_point, salsa
425
426    USE coriolis_mod,                                                          &
427        ONLY:  coriolis
428
429    USE cpulog,                                                                &
430        ONLY:  cpu_log, log_point, log_point_s
431
432    USE diffusion_s_mod,                                                       &
433        ONLY:  diffusion_s
434
435    USE diffusion_u_mod,                                                       &
436        ONLY:  diffusion_u
437
438    USE diffusion_v_mod,                                                       &
439        ONLY:  diffusion_v
440
441    USE diffusion_w_mod,                                                       &
442        ONLY:  diffusion_w
443
444    USE indices,                                                               &
445        ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
446               nzb, nzt, wall_flags_0
447
448    USE kinds
449
450    USE lsf_nudging_mod,                                                       &
451        ONLY:  ls_advec, nudge
452
453    USE module_interface,                                                      &
454        ONLY:  module_interface_actions, &
455               module_interface_non_advective_processes, &
456               module_interface_exchange_horiz, &
457               module_interface_prognostic_equations
458
459    USE ocean_mod,                                                             &
460        ONLY:  stokes_drift_terms, stokes_force,   &
461               wave_breaking, wave_breaking_term
462
463    USE plant_canopy_model_mod,                                                &
464        ONLY:  cthf, pcm_tendency
465
466#if defined( __rrtmg )
467    USE radiation_model_mod,                                                   &
468        ONLY:  radiation, radiation_tendency,                                  &
469               skip_time_do_radiation
470#endif
471
472    USE statistics,                                                            &
473        ONLY:  hom
474
475    USE subsidence_mod,                                                        &
476        ONLY:  subsidence
477
478    USE surface_mod,                                                           &
479        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
480                surf_usm_v
481
482    USE turbulence_closure_mod,                                                &
483        ONLY:  tcm_prognostic_equations
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 turbulence closure
1060          CALL tcm_prognostic_equations( i, j, i_omp_start, tn )
1061!
1062!--       Calculate prognostic equations for all other modules
1063          CALL module_interface_prognostic_equations( i, j, i_omp_start, tn )
1064
1065       ENDDO  ! loop over j
1066    ENDDO  ! loop over i
1067    !$OMP END PARALLEL
1068
1069
1070    CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1071
1072    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_cache', 'end' )
1073
1074 END SUBROUTINE prognostic_equations_cache
1075
1076
1077!------------------------------------------------------------------------------!
1078! Description:
1079! ------------
1080!> Version for vector machines
1081!------------------------------------------------------------------------------!
1082
1083 SUBROUTINE prognostic_equations_vector
1084
1085
1086    INTEGER(iwp) ::  i     !<
1087    INTEGER(iwp) ::  j     !<
1088    INTEGER(iwp) ::  k     !<
1089
1090    REAL(wp)     ::  sbt  !<
1091
1092
1093    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_vector', 'start' )
1094!
1095!-- Calculate non advective processes for all other modules
1096    CALL module_interface_non_advective_processes
1097!
1098!-- Module Inferface for exchange horiz after non_advective_processes but before
1099!-- advection. Therefore, non_advective_processes must not run for ghost points.
1100    CALL module_interface_exchange_horiz()
1101!
1102!-- u-velocity component
1103    CALL cpu_log( log_point(5), 'u-equation', 'start' )
1104
1105    !$ACC KERNELS PRESENT(tend)
1106    tend = 0.0_wp
1107    !$ACC END KERNELS
1108    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1109       IF ( ws_scheme_mom )  THEN
1110          CALL advec_u_ws
1111       ELSE
1112          CALL advec_u_pw
1113       ENDIF
1114    ELSE
1115       CALL advec_u_up
1116    ENDIF
1117    CALL diffusion_u
1118    CALL coriolis( 1 )
1119    IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
1120       CALL buoyancy( pt, 1 )
1121    ENDIF
1122
1123!
1124!-- Drag by plant canopy
1125    IF ( plant_canopy )  CALL pcm_tendency( 1 )
1126
1127!
1128!-- External pressure gradient
1129    IF ( dp_external )  THEN
1130       DO  i = nxlu, nxr
1131          DO  j = nys, nyn
1132             DO  k = dp_level_ind_b+1, nzt
1133                tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
1134             ENDDO
1135          ENDDO
1136       ENDDO
1137    ENDIF
1138
1139!
1140!-- Nudging
1141    IF ( nudging )  CALL nudge( simulated_time, 'u' )
1142
1143!
1144!-- Effect of Stokes drift (in ocean mode only)
1145    IF ( stokes_force )  CALL stokes_drift_terms( 1 )
1146
1147    CALL module_interface_actions( 'u-tendency' )
1148
1149!
1150!-- Prognostic equation for u-velocity component
1151    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1152    !$ACC PRESENT(u, tend, tu_m, u_init, rdf, wall_flags_0) &
1153    !$ACC PRESENT(tsc(2:5)) &
1154    !$ACC PRESENT(u_p)
1155    DO  i = nxlu, nxr
1156       DO  j = nys, nyn
1157          DO  k = nzb+1, nzt
1158             u_p(k,j,i) = u(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +          &
1159                                                 tsc(3) * tu_m(k,j,i) )          &
1160                                               - tsc(5) * rdf(k) *               &
1161                                                        ( u(k,j,i) - u_init(k) ) &
1162                                     ) * MERGE( 1.0_wp, 0.0_wp,                  &
1163                                                 BTEST( wall_flags_0(k,j,i), 1 ) &
1164                                              )
1165          ENDDO
1166       ENDDO
1167    ENDDO
1168
1169!
1170!-- Add turbulence generated by wave breaking (in ocean mode only)
1171    IF ( wave_breaking  .AND.                                                  &
1172         intermediate_timestep_count == intermediate_timestep_count_max )      &
1173    THEN
1174       CALL wave_breaking_term( 1 )
1175    ENDIF
1176
1177!
1178!-- Calculate tendencies for the next Runge-Kutta step
1179    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1180       IF ( intermediate_timestep_count == 1 )  THEN
1181          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1182          !$ACC PRESENT(tend, tu_m)
1183          DO  i = nxlu, nxr
1184             DO  j = nys, nyn
1185                DO  k = nzb+1, nzt
1186                   tu_m(k,j,i) = tend(k,j,i)
1187                ENDDO
1188             ENDDO
1189          ENDDO
1190       ELSEIF ( intermediate_timestep_count < &
1191                intermediate_timestep_count_max )  THEN
1192          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1193          !$ACC PRESENT(tend, tu_m)
1194          DO  i = nxlu, nxr
1195             DO  j = nys, nyn
1196                DO  k = nzb+1, nzt
1197                   tu_m(k,j,i) =    -9.5625_wp * tend(k,j,i)                   &
1198                                   + 5.3125_wp * tu_m(k,j,i)
1199                ENDDO
1200             ENDDO
1201          ENDDO
1202       ENDIF
1203    ENDIF
1204
1205    CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1206
1207!
1208!-- v-velocity component
1209    CALL cpu_log( log_point(6), 'v-equation', 'start' )
1210
1211    !$ACC KERNELS PRESENT(tend)
1212    tend = 0.0_wp
1213    !$ACC END KERNELS
1214    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1215       IF ( ws_scheme_mom )  THEN
1216          CALL advec_v_ws
1217       ELSE
1218          CALL advec_v_pw
1219       END IF
1220    ELSE
1221       CALL advec_v_up
1222    ENDIF
1223    CALL diffusion_v
1224    CALL coriolis( 2 )
1225
1226!
1227!-- Drag by plant canopy
1228    IF ( plant_canopy )  CALL pcm_tendency( 2 )
1229
1230!
1231!-- External pressure gradient
1232    IF ( dp_external )  THEN
1233       DO  i = nxl, nxr
1234          DO  j = nysv, nyn
1235             DO  k = dp_level_ind_b+1, nzt
1236                tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1237             ENDDO
1238          ENDDO
1239       ENDDO
1240    ENDIF
1241
1242!
1243!-- Nudging
1244    IF ( nudging )  CALL nudge( simulated_time, 'v' )
1245
1246!
1247!-- Effect of Stokes drift (in ocean mode only)
1248    IF ( stokes_force )  CALL stokes_drift_terms( 2 )
1249
1250    CALL module_interface_actions( 'v-tendency' )
1251
1252!
1253!-- Prognostic equation for v-velocity component
1254    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1255    !$ACC PRESENT(v, tend, tv_m, v_init, rdf, wall_flags_0) &
1256    !$ACC PRESENT(tsc(2:5)) &
1257    !$ACC PRESENT(v_p)
1258    DO  i = nxl, nxr
1259       DO  j = nysv, nyn
1260          DO  k = nzb+1, nzt
1261             v_p(k,j,i) = v(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1262                                                 tsc(3) * tv_m(k,j,i) )        &
1263                                               - tsc(5) * rdf(k) *             &
1264                                                      ( v(k,j,i) - v_init(k) ) &
1265                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1266                                                BTEST( wall_flags_0(k,j,i), 2 )&
1267                                              )
1268          ENDDO
1269       ENDDO
1270    ENDDO
1271
1272!
1273!-- Add turbulence generated by wave breaking (in ocean mode only)
1274    IF ( wave_breaking  .AND.                                                  &
1275         intermediate_timestep_count == intermediate_timestep_count_max )      &
1276    THEN
1277       CALL wave_breaking_term( 2 )
1278    ENDIF
1279
1280!
1281!-- Calculate tendencies for the next Runge-Kutta step
1282    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1283       IF ( intermediate_timestep_count == 1 )  THEN
1284          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1285          !$ACC PRESENT(tend, tv_m)
1286          DO  i = nxl, nxr
1287             DO  j = nysv, nyn
1288                DO  k = nzb+1, nzt
1289                   tv_m(k,j,i) = tend(k,j,i)
1290                ENDDO
1291             ENDDO
1292          ENDDO
1293       ELSEIF ( intermediate_timestep_count < &
1294                intermediate_timestep_count_max )  THEN
1295          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1296          !$ACC PRESENT(tend, tv_m)
1297          DO  i = nxl, nxr
1298             DO  j = nysv, nyn
1299                DO  k = nzb+1, nzt
1300                   tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1301                                  + 5.3125_wp * tv_m(k,j,i)
1302                ENDDO
1303             ENDDO
1304          ENDDO
1305       ENDIF
1306    ENDIF
1307
1308    CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1309
1310!
1311!-- w-velocity component
1312    CALL cpu_log( log_point(7), 'w-equation', 'start' )
1313
1314    !$ACC KERNELS PRESENT(tend)
1315    tend = 0.0_wp
1316    !$ACC END KERNELS
1317    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1318       IF ( ws_scheme_mom )  THEN
1319          CALL advec_w_ws
1320       ELSE
1321          CALL advec_w_pw
1322       ENDIF
1323    ELSE
1324       CALL advec_w_up
1325    ENDIF
1326    CALL diffusion_w
1327    CALL coriolis( 3 )
1328
1329    IF ( .NOT. neutral )  THEN
1330       IF ( ocean_mode )  THEN
1331          CALL buoyancy( rho_ocean, 3 )
1332       ELSE
1333          IF ( .NOT. humidity )  THEN
1334             CALL buoyancy( pt, 3 )
1335          ELSE
1336             CALL buoyancy( vpt, 3 )
1337          ENDIF
1338       ENDIF
1339    ENDIF
1340
1341!
1342!-- Drag by plant canopy
1343    IF ( plant_canopy )  CALL pcm_tendency( 3 )
1344
1345!
1346!-- Effect of Stokes drift (in ocean mode only)
1347    IF ( stokes_force )  CALL stokes_drift_terms( 3 )
1348
1349    CALL module_interface_actions( 'w-tendency' )
1350
1351!
1352!-- Prognostic equation for w-velocity component
1353    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1354    !$ACC PRESENT(w, tend, tw_m, v_init, rdf, wall_flags_0) &
1355    !$ACC PRESENT(tsc(2:5)) &
1356    !$ACC PRESENT(w_p)
1357    DO  i = nxl, nxr
1358       DO  j = nys, nyn
1359          DO  k = nzb+1, nzt-1
1360             w_p(k,j,i) = w(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1361                                                 tsc(3) * tw_m(k,j,i) )        &
1362                                               - tsc(5) * rdf(k) * w(k,j,i)    &
1363                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1364                                                BTEST( wall_flags_0(k,j,i), 3 )&
1365                                              )
1366          ENDDO
1367       ENDDO
1368    ENDDO
1369
1370!
1371!-- Calculate tendencies for the next Runge-Kutta step
1372    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1373       IF ( intermediate_timestep_count == 1 )  THEN
1374          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1375          !$ACC PRESENT(tend, tw_m)
1376          DO  i = nxl, nxr
1377             DO  j = nys, nyn
1378                DO  k = nzb+1, nzt-1
1379                   tw_m(k,j,i) = tend(k,j,i)
1380                ENDDO
1381             ENDDO
1382          ENDDO
1383       ELSEIF ( intermediate_timestep_count < &
1384                intermediate_timestep_count_max )  THEN
1385          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1386          !$ACC PRESENT(tend, tw_m)
1387          DO  i = nxl, nxr
1388             DO  j = nys, nyn
1389                DO  k = nzb+1, nzt-1
1390                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1391                                  + 5.3125_wp * tw_m(k,j,i)
1392                ENDDO
1393             ENDDO
1394          ENDDO
1395       ENDIF
1396    ENDIF
1397
1398    CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1399
1400
1401!
1402!-- If required, compute prognostic equation for potential temperature
1403    IF ( .NOT. neutral )  THEN
1404
1405       CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1406
1407!
1408!--    pt-tendency terms with communication
1409       sbt = tsc(2)
1410       IF ( scalar_advec == 'bc-scheme' )  THEN
1411
1412          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1413!
1414!--          Bott-Chlond scheme always uses Euler time step. Thus:
1415             sbt = 1.0_wp
1416          ENDIF
1417          tend = 0.0_wp
1418          CALL advec_s_bc( pt, 'pt' )
1419
1420       ENDIF
1421
1422!
1423!--    pt-tendency terms with no communication
1424       IF ( scalar_advec /= 'bc-scheme' )  THEN
1425          !$ACC KERNELS PRESENT(tend)
1426          tend = 0.0_wp
1427          !$ACC END KERNELS
1428          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1429             IF ( ws_scheme_sca )  THEN
1430                CALL advec_s_ws( pt, 'pt' )
1431             ELSE
1432                CALL advec_s_pw( pt )
1433             ENDIF
1434          ELSE
1435             CALL advec_s_up( pt )
1436          ENDIF
1437       ENDIF
1438
1439       CALL diffusion_s( pt,                                                   &
1440                         surf_def_h(0)%shf, surf_def_h(1)%shf,                 &
1441                         surf_def_h(2)%shf,                                    &
1442                         surf_lsm_h%shf,    surf_usm_h%shf,                    &
1443                         surf_def_v(0)%shf, surf_def_v(1)%shf,                 &
1444                         surf_def_v(2)%shf, surf_def_v(3)%shf,                 &
1445                         surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,                 &
1446                         surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,                 &
1447                         surf_usm_v(0)%shf, surf_usm_v(1)%shf,                 &
1448                         surf_usm_v(2)%shf, surf_usm_v(3)%shf )
1449
1450!
1451!--    Consideration of heat sources within the plant canopy
1452       IF ( plant_canopy  .AND.                                          &
1453            (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
1454          CALL pcm_tendency( 4 )
1455       ENDIF
1456
1457!
1458!--    Large scale advection
1459       IF ( large_scale_forcing )  THEN
1460          CALL ls_advec( simulated_time, 'pt' )
1461       ENDIF
1462
1463!
1464!--    Nudging
1465       IF ( nudging )  CALL nudge( simulated_time, 'pt' )
1466
1467!
1468!--    If required compute influence of large-scale subsidence/ascent
1469       IF ( large_scale_subsidence  .AND.                                      &
1470            .NOT. use_subsidence_tendencies )  THEN
1471          CALL subsidence( tend, pt, pt_init, 2 )
1472       ENDIF
1473
1474#if defined( __rrtmg )
1475!
1476!--    If required, add tendency due to radiative heating/cooling
1477       IF ( radiation  .AND.                                                   &
1478            simulated_time > skip_time_do_radiation )  THEN
1479            CALL radiation_tendency ( tend )
1480       ENDIF
1481#endif
1482
1483       CALL module_interface_actions( 'pt-tendency' )
1484
1485!
1486!--    Prognostic equation for potential temperature
1487       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1488       !$ACC PRESENT(pt, tend, tpt_m, wall_flags_0) &
1489       !$ACC PRESENT(pt_init, rdf_sc, ptdf_x, ptdf_y) &
1490       !$ACC PRESENT(tsc(3:5)) &
1491       !$ACC PRESENT(pt_p)
1492       DO  i = nxl, nxr
1493          DO  j = nys, nyn
1494             DO  k = nzb+1, nzt
1495                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +      &
1496                                                   tsc(3) * tpt_m(k,j,i) )     &
1497                                                 - tsc(5) *                    &
1498                                                   ( pt(k,j,i) - pt_init(k) ) *&
1499                                          ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )&
1500                                          )                                    &
1501                                   * MERGE( 1.0_wp, 0.0_wp,                    &
1502                                             BTEST( wall_flags_0(k,j,i), 0 )   &
1503                                          )
1504             ENDDO
1505          ENDDO
1506       ENDDO
1507!
1508!--    Calculate tendencies for the next Runge-Kutta step
1509       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1510          IF ( intermediate_timestep_count == 1 )  THEN
1511             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1512             !$ACC PRESENT(tend, tpt_m)
1513             DO  i = nxl, nxr
1514                DO  j = nys, nyn
1515                   DO  k = nzb+1, nzt
1516                      tpt_m(k,j,i) = tend(k,j,i)
1517                   ENDDO
1518                ENDDO
1519             ENDDO
1520          ELSEIF ( intermediate_timestep_count < &
1521                   intermediate_timestep_count_max )  THEN
1522             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1523             !$ACC PRESENT(tend, tpt_m)
1524             DO  i = nxl, nxr
1525                DO  j = nys, nyn
1526                   DO  k = nzb+1, nzt
1527                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
1528                                        5.3125_wp * tpt_m(k,j,i)
1529                   ENDDO
1530                ENDDO
1531             ENDDO
1532          ENDIF
1533       ENDIF
1534
1535       CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1536
1537    ENDIF
1538
1539!
1540!-- If required, compute prognostic equation for total water content
1541    IF ( humidity )  THEN
1542
1543       CALL cpu_log( log_point(29), 'q-equation', 'start' )
1544
1545!
1546!--    Scalar/q-tendency terms with communication
1547       sbt = tsc(2)
1548       IF ( scalar_advec == 'bc-scheme' )  THEN
1549
1550          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1551!
1552!--          Bott-Chlond scheme always uses Euler time step. Thus:
1553             sbt = 1.0_wp
1554          ENDIF
1555          tend = 0.0_wp
1556          CALL advec_s_bc( q, 'q' )
1557
1558       ENDIF
1559
1560!
1561!--    Scalar/q-tendency terms with no communication
1562       IF ( scalar_advec /= 'bc-scheme' )  THEN
1563          tend = 0.0_wp
1564          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1565             IF ( ws_scheme_sca )  THEN
1566                CALL advec_s_ws( q, 'q' )
1567             ELSE
1568                CALL advec_s_pw( q )
1569             ENDIF
1570          ELSE
1571             CALL advec_s_up( q )
1572          ENDIF
1573       ENDIF
1574
1575       CALL diffusion_s( q,                                                    &
1576                         surf_def_h(0)%qsws, surf_def_h(1)%qsws,               &
1577                         surf_def_h(2)%qsws,                                   &
1578                         surf_lsm_h%qsws,    surf_usm_h%qsws,                  &
1579                         surf_def_v(0)%qsws, surf_def_v(1)%qsws,               &
1580                         surf_def_v(2)%qsws, surf_def_v(3)%qsws,               &
1581                         surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,               &
1582                         surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,               &
1583                         surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,               &
1584                         surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
1585
1586!
1587!--    Sink or source of humidity due to canopy elements
1588       IF ( plant_canopy ) CALL pcm_tendency( 5 )
1589
1590!
1591!--    Large scale advection
1592       IF ( large_scale_forcing )  THEN
1593          CALL ls_advec( simulated_time, 'q' )
1594       ENDIF
1595
1596!
1597!--    Nudging
1598       IF ( nudging )  CALL nudge( simulated_time, 'q' )
1599
1600!
1601!--    If required compute influence of large-scale subsidence/ascent
1602       IF ( large_scale_subsidence  .AND.                                      &
1603            .NOT. use_subsidence_tendencies )  THEN
1604         CALL subsidence( tend, q, q_init, 3 )
1605       ENDIF
1606
1607       CALL module_interface_actions( 'q-tendency' )
1608
1609!
1610!--    Prognostic equation for total water content
1611       DO  i = nxl, nxr
1612          DO  j = nys, nyn
1613             DO  k = nzb+1, nzt
1614                q_p(k,j,i) = q(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
1615                                                 tsc(3) * tq_m(k,j,i) )        &
1616                                               - tsc(5) * rdf_sc(k) *          &
1617                                                      ( q(k,j,i) - q_init(k) ) &
1618                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
1619                                                   BTEST( wall_flags_0(k,j,i), 0 ) &
1620                                                 )
1621                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
1622             ENDDO
1623          ENDDO
1624       ENDDO
1625
1626!
1627!--    Calculate tendencies for the next Runge-Kutta step
1628       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1629          IF ( intermediate_timestep_count == 1 )  THEN
1630             DO  i = nxl, nxr
1631                DO  j = nys, nyn
1632                   DO  k = nzb+1, nzt
1633                      tq_m(k,j,i) = tend(k,j,i)
1634                   ENDDO
1635                ENDDO
1636             ENDDO
1637          ELSEIF ( intermediate_timestep_count < &
1638                   intermediate_timestep_count_max )  THEN
1639             DO  i = nxl, nxr
1640                DO  j = nys, nyn
1641                   DO  k = nzb+1, nzt
1642                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
1643                                     + 5.3125_wp * tq_m(k,j,i)
1644                   ENDDO
1645                ENDDO
1646             ENDDO
1647          ENDIF
1648       ENDIF
1649
1650       CALL cpu_log( log_point(29), 'q-equation', 'stop' )
1651
1652    ENDIF
1653!
1654!-- If required, compute prognostic equation for scalar
1655    IF ( passive_scalar )  THEN
1656
1657       CALL cpu_log( log_point(66), 's-equation', 'start' )
1658
1659!
1660!--    Scalar/q-tendency terms with communication
1661       sbt = tsc(2)
1662       IF ( scalar_advec == 'bc-scheme' )  THEN
1663
1664          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1665!
1666!--          Bott-Chlond scheme always uses Euler time step. Thus:
1667             sbt = 1.0_wp
1668          ENDIF
1669          tend = 0.0_wp
1670          CALL advec_s_bc( s, 's' )
1671
1672       ENDIF
1673
1674!
1675!--    Scalar/q-tendency terms with no communication
1676       IF ( scalar_advec /= 'bc-scheme' )  THEN
1677          tend = 0.0_wp
1678          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1679             IF ( ws_scheme_sca )  THEN
1680                CALL advec_s_ws( s, 's' )
1681             ELSE
1682                CALL advec_s_pw( s )
1683             ENDIF
1684          ELSE
1685             CALL advec_s_up( s )
1686          ENDIF
1687       ENDIF
1688
1689       CALL diffusion_s( s,                                                    &
1690                         surf_def_h(0)%ssws, surf_def_h(1)%ssws,               &
1691                         surf_def_h(2)%ssws,                                   &
1692                         surf_lsm_h%ssws,    surf_usm_h%ssws,                  &
1693                         surf_def_v(0)%ssws, surf_def_v(1)%ssws,               &
1694                         surf_def_v(2)%ssws, surf_def_v(3)%ssws,               &
1695                         surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,               &
1696                         surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,               &
1697                         surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,               &
1698                         surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1699
1700!
1701!--    Sink or source of humidity due to canopy elements
1702       IF ( plant_canopy ) CALL pcm_tendency( 7 )
1703
1704!
1705!--    Large scale advection. Not implemented for scalars so far.
1706!        IF ( large_scale_forcing )  THEN
1707!           CALL ls_advec( simulated_time, 'q' )
1708!        ENDIF
1709
1710!
1711!--    Nudging. Not implemented for scalars so far.
1712!        IF ( nudging )  CALL nudge( simulated_time, 'q' )
1713
1714!
1715!--    If required compute influence of large-scale subsidence/ascent.
1716!--    Not implemented for scalars so far.
1717       IF ( large_scale_subsidence  .AND.                                      &
1718            .NOT. use_subsidence_tendencies  .AND.                             &
1719            .NOT. large_scale_forcing )  THEN
1720         CALL subsidence( tend, s, s_init, 3 )
1721       ENDIF
1722
1723       CALL module_interface_actions( 's-tendency' )
1724
1725!
1726!--    Prognostic equation for total water content
1727       DO  i = nxl, nxr
1728          DO  j = nys, nyn
1729             DO  k = nzb+1, nzt
1730                s_p(k,j,i) = s(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
1731                                                 tsc(3) * ts_m(k,j,i) )        &
1732                                               - tsc(5) * rdf_sc(k) *          &
1733                                                    ( s(k,j,i) - s_init(k) )   &
1734                                        )                                      &
1735                                   * MERGE( 1.0_wp, 0.0_wp,                    &
1736                                             BTEST( wall_flags_0(k,j,i), 0 )   &
1737                                          )
1738                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
1739             ENDDO
1740          ENDDO
1741       ENDDO
1742
1743!
1744!--    Calculate tendencies for the next Runge-Kutta step
1745       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1746          IF ( intermediate_timestep_count == 1 )  THEN
1747             DO  i = nxl, nxr
1748                DO  j = nys, nyn
1749                   DO  k = nzb+1, nzt
1750                      ts_m(k,j,i) = tend(k,j,i)
1751                   ENDDO
1752                ENDDO
1753             ENDDO
1754          ELSEIF ( intermediate_timestep_count < &
1755                   intermediate_timestep_count_max )  THEN
1756             DO  i = nxl, nxr
1757                DO  j = nys, nyn
1758                   DO  k = nzb+1, nzt
1759                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
1760                                     + 5.3125_wp * ts_m(k,j,i)
1761                   ENDDO
1762                ENDDO
1763             ENDDO
1764          ENDIF
1765       ENDIF
1766
1767       CALL cpu_log( log_point(66), 's-equation', 'stop' )
1768
1769    ENDIF
1770
1771!
1772!-- Calculate prognostic equations for turbulence closure
1773    CALL tcm_prognostic_equations()
1774!
1775!-- Calculate prognostic equations for all other modules
1776    CALL module_interface_prognostic_equations()
1777
1778    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_vector', 'end' )
1779
1780 END SUBROUTINE prognostic_equations_vector
1781
1782
1783 END MODULE prognostic_equations_mod
Note: See TracBrowser for help on using the repository browser.