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

Last change on this file since 3961 was 3956, checked in by monakurppa, 5 years ago

Remove salsa calls from prognostic_equations and correct a bug in the salsa deposition for urban and land surface models

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