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

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

Moved chem_prognostic_equations into module_interface

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