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

Last change on this file since 3885 was 3885, checked in by kanani, 6 years ago

restructure/add location/debug messages

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