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

Last change on this file since 3881 was 3881, checked in by suehring, 5 years ago

Make level 3 initialization of urban-surfaces consistent to input data standard; revise flagging of ground-floor level surfaces at buidlings; bugfixes in level 3 initialization of albedo; bugfix in OpenMP directive; check for zero output timestep in surface output; assure that Obukhov lenght does not become zero

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