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

Last change on this file since 3822 was 3820, checked in by forkel, 6 years ago

renaming of get_mechanismname, do_emiss and do_depo, sorting in chem_modules

  • Property svn:keywords set to Id
File size: 104.8 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 3820 2019-03-27 11:53:41Z hellstea $
27! renamed do_depo to deposition_dry (ecc)
28!
29! 3797 2019-03-15 11:15:38Z forkel
30! Call chem_integegrate in OpenMP loop   (ketelsen)
31!
32!
33! 3771 2019-02-28 12:19:33Z raasch
34! preprocessor directivs fro rrtmg added
35!
36! 3761 2019-02-25 15:31:42Z raasch
37! unused variable removed
38!
39! 3719 2019-02-06 13:10:18Z kanani
40! Cleaned up chemistry cpu measurements
41!
42! 3684 2019-01-20 20:20:58Z knoop
43! OpenACC port for SPEC
44!
45! 3589 2018-11-30 15:09:51Z suehring
46! Move the control parameter "salsa" from salsa_mod to control_parameters
47! (M. Kurppa)
48!
49! 3582 2018-11-29 19:16:36Z suehring
50! Implementation of a new aerosol module salsa.
51! Remove cpu-logs from i,j loop in cache version.
52!
53! 3458 2018-10-30 14:51:23Z kanani
54! remove duplicate USE chem_modules
55! from chemistry branch r3443, banzhafs, basit:
56! chem_depo call introduced
57! code added for decycling chemistry
58!
59! 3386 2018-10-19 16:28:22Z gronemeier
60! Renamed tcm_prognostic to tcm_prognostic_equations
61!
62! 3355 2018-10-16 14:03:34Z knoop
63! (from branch resler)
64! Fix for chemistry call
65!
66! 3302 2018-10-03 02:39:40Z raasch
67! Stokes drift + wave breaking term added
68!
69! 3298 2018-10-02 12:21:11Z kanani
70! Code added for decycling chemistry (basit)
71!
72! 3294 2018-10-01 02:37:10Z raasch
73! changes concerning modularization of ocean option
74!
75! 3274 2018-09-24 15:42:55Z knoop
76! Modularization of all bulk cloud physics code components
77!
78! 3241 2018-09-12 15:02:00Z raasch
79! omp_get_thread_num now declared in omp directive
80!
81! 3183 2018-07-27 14:25:55Z suehring
82! Remove unused variables from USE statements
83!
84! 3182 2018-07-27 13:36:03Z suehring
85! Revise recent bugfix for nesting
86!
87! 3021 2018-05-16 08:14:20Z maronga
88! Bugfix in IF clause for nesting
89!
90! 3014 2018-05-09 08:42:38Z maronga
91! Fixed a bug in the IF condition to call pcm_tendency in case of
92! potential temperature
93!
94! 2815 2018-02-19 11:29:57Z kanani
95! Rename chem_tendency to chem_prognostic_equations,
96! implement vector version for air chemistry
97!
98! 2766 2018-01-22 17:17:47Z kanani
99! Removed preprocessor directive __chem
100!
101! 2746 2018-01-15 12:06:04Z suehring
102! Move flag plant canopy to modules
103!
104! 2719 2018-01-02 09:02:06Z maronga
105! Bugfix for last change.
106!
107! 2718 2018-01-02 08:49:38Z maronga
108! Corrected "Former revisions" section
109!
110! 2696 2017-12-14 17:12:51Z kanani
111! - Change in file header (GPL part)
112! - Moved TKE equation to tcm_prognostic (TG)
113! - Added switch for chemical reactions (RF, FK)
114! - Implementation of chemistry module (RF, BK, FK)
115!
116! 2563 2017-10-19 15:36:10Z Giersch
117! Variable wind_turbine moved to module control_parameters
118!
119! 2320 2017-07-21 12:47:43Z suehring
120! Modularize large-scale forcing and nudging
121!
122! 2292 2017-06-20 09:51:42Z schwenkel
123! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
124! includes two more prognostic equations for cloud drop concentration (nc) 
125! and cloud water content (qc).
126!
127! 2261 2017-06-08 14:25:57Z raasch
128! bugfix for r2232: openmp directives removed
129!
130! 2233 2017-05-30 18:08:54Z suehring
131!
132! 2232 2017-05-30 17:47:52Z suehring
133! Adjutst to new surface-type structure. Remove call for usm_wall_heat_flux,
134! which is realized directly in diffusion_s now.
135!
136! 2192 2017-03-22 04:14:10Z raasch
137! Bugfix for misplaced and missing openMP directives from r2155
138!
139! 2155 2017-02-21 09:57:40Z hoffmann
140! Bugfix in the calculation of microphysical quantities on ghost points.
141!
142! 2118 2017-01-17 16:38:49Z raasch
143! OpenACC version of subroutine removed
144!
145! 2031 2016-10-21 15:11:58Z knoop
146! renamed variable rho to rho_ocean
147!
148! 2011 2016-09-19 17:29:57Z kanani
149! Flag urban_surface is now defined in module control_parameters.
150!
151! 2007 2016-08-24 15:47:17Z kanani
152! Added pt tendency calculation based on energy balance at urban surfaces
153! (new urban surface model)
154!
155! 2000 2016-08-20 18:09:15Z knoop
156! Forced header and separation lines into 80 columns
157!
158! 1976 2016-07-27 13:28:04Z maronga
159! Simplied calls to radiation model
160!
161! 1960 2016-07-12 16:34:24Z suehring
162! Separate humidity and passive scalar
163!
164! 1914 2016-05-26 14:44:07Z witha
165! Added calls for wind turbine model
166!
167! 1873 2016-04-18 14:50:06Z maronga
168! Module renamed (removed _mod)
169!
170! 1850 2016-04-08 13:29:27Z maronga
171! Module renamed
172!
173! 1826 2016-04-07 12:01:39Z maronga
174! Renamed canopy model calls.
175!
176! 1822 2016-04-07 07:49:42Z hoffmann
177! Kessler microphysics scheme moved to microphysics.
178!
179! 1757 2016-02-22 15:49:32Z maronga
180!
181! 1691 2015-10-26 16:17:44Z maronga
182! Added optional model spin-up without radiation / land surface model calls.
183! Formatting corrections.
184!
185! 1682 2015-10-07 23:56:08Z knoop
186! Code annotations made doxygen readable
187!
188! 1585 2015-04-30 07:05:52Z maronga
189! Added call for temperature tendency calculation due to radiative flux divergence
190!
191! 1517 2015-01-07 19:12:25Z hoffmann
192! advec_s_bc_mod addded, since advec_s_bc is now a module
193!
194! 1484 2014-10-21 10:53:05Z kanani
195! Changes due to new module structure of the plant canopy model:
196! parameters cthf and plant_canopy moved to module plant_canopy_model_mod.
197! Removed double-listing of use_upstream_for_tke in ONLY-list of module
198! control_parameters
199!
200! 1409 2014-05-23 12:11:32Z suehring
201! Bugfix: i_omp_start changed for advec_u_ws at left inflow and outflow boundary.
202! This ensures that left-hand side fluxes are also calculated for nxl in that
203! case, even though the solution at nxl is overwritten in boundary_conds()
204!
205! 1398 2014-05-07 11:15:00Z heinze
206! Rayleigh-damping for horizontal velocity components changed: instead of damping
207! against ug and vg, damping against u_init and v_init is used to allow for a
208! homogenized treatment in case of nudging
209!
210! 1380 2014-04-28 12:40:45Z heinze
211! Change order of calls for scalar prognostic quantities:
212! ls_advec -> nudging -> subsidence since initial profiles
213!
214! 1374 2014-04-25 12:55:07Z raasch
215! missing variables added to ONLY lists
216!
217! 1365 2014-04-22 15:03:56Z boeske
218! Calls of ls_advec for large scale advection added,
219! subroutine subsidence is only called if use_subsidence_tendencies = .F.,
220! new argument ls_index added to the calls of subsidence
221! +ls_index
222!
223! 1361 2014-04-16 15:17:48Z hoffmann
224! Two-moment microphysics moved to the start of prognostic equations. This makes
225! the 3d arrays for tend_q, tend_qr, tend_pt and tend_pt redundant.
226! Additionally, it is allowed to call the microphysics just once during the time
227! step (not at each sub-time step).
228!
229! Two-moment cloud physics added for vector and accelerator optimization.
230!
231! 1353 2014-04-08 15:21:23Z heinze
232! REAL constants provided with KIND-attribute
233!
234! 1337 2014-03-25 15:11:48Z heinze
235! Bugfix: REAL constants provided with KIND-attribute
236!
237! 1332 2014-03-25 11:59:43Z suehring
238! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke
239!
240! 1330 2014-03-24 17:29:32Z suehring
241! In case of SGS-particle velocity advection of TKE is also allowed with
242! dissipative 5th-order scheme.
243!
244! 1320 2014-03-20 08:40:49Z raasch
245! ONLY-attribute added to USE-statements,
246! kind-parameters added to all INTEGER and REAL declaration statements,
247! kinds are defined in new module kinds,
248! old module precision_kind is removed,
249! revision history before 2012 removed,
250! comment fields (!:) to be used for variable explanations added to
251! all variable declaration statements
252!
253! 1318 2014-03-17 13:35:16Z raasch
254! module interfaces removed
255!
256! 1257 2013-11-08 15:18:40Z raasch
257! openacc loop vector clauses removed, independent clauses added
258!
259! 1246 2013-11-01 08:59:45Z heinze
260! enable nudging also for accelerator version
261!
262! 1241 2013-10-30 11:36:58Z heinze
263! usage of nudging enabled (so far not implemented for accelerator version)
264!
265! 1179 2013-06-14 05:57:58Z raasch
266! two arguments removed from routine buoyancy, ref_state updated on device
267!
268! 1128 2013-04-12 06:19:32Z raasch
269! those parts requiring global communication moved to time_integration,
270! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
271! j_north
272!
273! 1115 2013-03-26 18:16:16Z hoffmann
274! optimized cloud physics: calculation of microphysical tendencies transfered
275! to microphysics.f90; qr and nr are only calculated if precipitation is required
276!
277! 1111 2013-03-08 23:54:10Z raasch
278! update directives for prognostic quantities removed
279!
280! 1106 2013-03-04 05:31:38Z raasch
281! small changes in code formatting
282!
283! 1092 2013-02-02 11:24:22Z raasch
284! unused variables removed
285!
286! 1053 2012-11-13 17:11:03Z hoffmann
287! implementation of two new prognostic equations for rain drop concentration (nr)
288! and rain water content (qr)
289!
290! currently, only available for cache loop optimization
291!
292! 1036 2012-10-22 13:43:42Z raasch
293! code put under GPL (PALM 3.9)
294!
295! 1019 2012-09-28 06:46:45Z raasch
296! non-optimized version of prognostic_equations removed
297!
298! 1015 2012-09-27 09:23:24Z raasch
299! new branch prognostic_equations_acc
300! OpenACC statements added + code changes required for GPU optimization
301!
302! 1001 2012-09-13 14:08:46Z raasch
303! all actions concerning leapfrog- and upstream-spline-scheme removed
304!
305! 978 2012-08-09 08:28:32Z fricke
306! km_damp_x and km_damp_y removed in calls of diffusion_u and diffusion_v
307! add ptdf_x, ptdf_y for damping the potential temperature at the inflow
308! boundary in case of non-cyclic lateral boundaries
309! Bugfix: first thread index changes for WS-scheme at the inflow
310!
311! 940 2012-07-09 14:31:00Z raasch
312! temperature equation can be switched off
313!
314! Revision 1.1  2000/04/13 14:56:27  schroeter
315! Initial revision
316!
317!
318! Description:
319! ------------
320!> Solving the prognostic equations.
321!------------------------------------------------------------------------------!
322 MODULE prognostic_equations_mod
323
324
325    USE advec_s_bc_mod,                                                        &
326        ONLY:  advec_s_bc
327
328    USE advec_s_pw_mod,                                                        &
329        ONLY:  advec_s_pw
330
331    USE advec_s_up_mod,                                                        &
332        ONLY:  advec_s_up
333
334    USE advec_u_pw_mod,                                                        &
335        ONLY:  advec_u_pw
336
337    USE advec_u_up_mod,                                                        &
338        ONLY:  advec_u_up
339
340    USE advec_v_pw_mod,                                                        &
341        ONLY:  advec_v_pw
342
343    USE advec_v_up_mod,                                                        &
344        ONLY:  advec_v_up
345
346    USE advec_w_pw_mod,                                                        &
347        ONLY:  advec_w_pw
348
349    USE advec_w_up_mod,                                                        &
350        ONLY:  advec_w_up
351
352    USE advec_ws,                                                              &
353        ONLY:  advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws
354
355    USE arrays_3d,                                                             &
356        ONLY:  diss_l_e, diss_l_nc, diss_l_nr, diss_l_pt, diss_l_q, diss_l_qc, &
357               diss_l_qr, diss_l_s, diss_l_sa, diss_s_e, diss_s_nc, diss_s_nr, &
358               diss_s_pt, diss_s_q, diss_s_qc, diss_s_qr, diss_s_s, diss_s_sa, &
359               e, e_p, flux_s_e, flux_s_nc, flux_s_nr, flux_s_pt, flux_s_q,    &
360               flux_s_qc, flux_s_qr, flux_s_s, flux_s_sa, flux_l_e, flux_l_nc, &
361               flux_l_nr, flux_l_pt, flux_l_q, flux_l_qc, flux_l_qr, flux_l_s, &
362               flux_l_sa, nc, nc_p, nr, nr_p, pt, ptdf_x, ptdf_y, pt_init,     &
363               pt_p, prho, q, q_init, q_p, qc, qc_p, qr, qr_p, rdf, rdf_sc,    &
364               ref_state, rho_ocean, s, s_init, s_p, tend, te_m, tnc_m,        &
365               tnr_m, tpt_m, tq_m, tqc_m, tqr_m, ts_m, tu_m, tv_m, tw_m, u,    &
366               ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
367
368    USE bulk_cloud_model_mod,                                                  &
369        ONLY:  call_microphysics_at_all_substeps, bulk_cloud_model,            &
370               bcm_actions, microphysics_sat_adjust,                           &
371               microphysics_morrison, microphysics_seifert
372
373    USE buoyancy_mod,                                                          &
374        ONLY:  buoyancy
375           
376    USE chem_modules,                                                          &
377        ONLY:  call_chem_at_all_substeps, chem_gasphase_on, cs_name,           &
378               deposition_dry
379
380    USE chem_photolysis_mod,                                                   &
381        ONLY:  photolysis_control
382
383    USE chemistry_model_mod,                                                   &
384        ONLY:  chem_boundary_conds, chem_depo, chem_integrate,                 &
385               chem_prognostic_equations, chem_species,                        &
386               nspec, nvar, spc_names
387
388    USE control_parameters,                                                    &
389        ONLY:  air_chemistry, constant_diffusion,                              &
390               dp_external, dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d,    &
391               humidity, intermediate_timestep_count,                          &
392               intermediate_timestep_count_max, large_scale_forcing,           &
393               large_scale_subsidence, neutral, nudging,                       &
394               ocean_mode, passive_scalar, plant_canopy, pt_reference,         &
395               scalar_advec, scalar_advec, simulated_time, sloping_surface,    &
396               timestep_scheme, tsc, use_subsidence_tendencies,                &
397               use_upstream_for_tke, wind_turbine, ws_scheme_mom,              & 
398               ws_scheme_sca, urban_surface, land_surface,                     &
399               time_since_reference_point, salsa
400
401    USE coriolis_mod,                                                          &
402        ONLY:  coriolis
403
404    USE cpulog,                                                                &
405        ONLY:  cpu_log, log_point, log_point_s
406
407    USE diffusion_s_mod,                                                       &
408        ONLY:  diffusion_s
409
410    USE diffusion_u_mod,                                                       &
411        ONLY:  diffusion_u
412
413    USE diffusion_v_mod,                                                       &
414        ONLY:  diffusion_v
415
416    USE diffusion_w_mod,                                                       &
417        ONLY:  diffusion_w
418
419    USE indices,                                                               &
420        ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
421               nzb, nzt, wall_flags_0
422
423    USE kinds
424
425    USE lsf_nudging_mod,                                                       &
426        ONLY:  ls_advec, nudge
427
428    USE module_interface,                                                      &
429        ONLY:  module_interface_actions
430
431    USE ocean_mod,                                                             &
432        ONLY:  ocean_prognostic_equations, stokes_drift_terms, stokes_force,   &
433               wave_breaking, wave_breaking_term
434
435    USE plant_canopy_model_mod,                                                &
436        ONLY:  cthf, pcm_tendency
437
438#if defined( __rrtmg )
439    USE radiation_model_mod,                                                   &
440        ONLY:  radiation, radiation_tendency,                                  &
441               skip_time_do_radiation
442#endif
443                         
444    USE salsa_mod,                                                             &
445        ONLY:  aerosol_mass, aerosol_number, dt_salsa, last_salsa_time, nbins, &
446               ncc_tot, ngast, salsa_boundary_conds, salsa_diagnostics,        &
447               salsa_driver, salsa_gas, salsa_gases_from_chem, salsa_tendency, &
448               skip_time_do_salsa
449               
450    USE salsa_util_mod,                                                        &
451        ONLY:  sums_salsa_ws_l                             
452   
453    USE statistics,                                                            &
454        ONLY:  hom
455
456    USE subsidence_mod,                                                        &
457        ONLY:  subsidence
458
459    USE surface_mod,                                                           &
460        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
461                surf_usm_v
462
463    USE turbulence_closure_mod,                                                &
464        ONLY:  tcm_prognostic_equations
465
466    USE wind_turbine_model_mod,                                                &
467        ONLY:  wtm_tendencies
468
469
470    PRIVATE
471    PUBLIC prognostic_equations_cache, prognostic_equations_vector
472
473    INTERFACE prognostic_equations_cache
474       MODULE PROCEDURE prognostic_equations_cache
475    END INTERFACE prognostic_equations_cache
476
477    INTERFACE prognostic_equations_vector
478       MODULE PROCEDURE prognostic_equations_vector
479    END INTERFACE prognostic_equations_vector
480
481
482 CONTAINS
483
484
485!------------------------------------------------------------------------------!
486! Description:
487! ------------
488!> Version with one optimized loop over all equations. It is only allowed to
489!> be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
490!>
491!> Here the calls of most subroutines are embedded in two DO loops over i and j,
492!> so communication between CPUs is not allowed (does not make sense) within
493!> these loops.
494!>
495!> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
496!------------------------------------------------------------------------------!
497
498 SUBROUTINE prognostic_equations_cache
499
500
501    IMPLICIT NONE
502   
503    INTEGER(iwp) ::  b                   !< index for aerosol size bins (salsa)
504    INTEGER(iwp) ::  c                   !< index for chemical compounds (salsa)
505    INTEGER(iwp) ::  g                   !< index for gaseous compounds (salsa)
506    INTEGER(iwp) ::  i                   !<
507    INTEGER(iwp) ::  i_omp_start         !<
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    INTEGER(iwp) ::  lsp
515    INTEGER(iwp) ::  lsp_usr             !< lsp running index for chem spcs
516
517
518!
519!-- Time measurement can only be performed for the whole set of equations
520    CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
521
522!
523!-- Calculation of chemical reactions. This is done outside of main loop,
524!-- since exchange of ghost points is required after this update of the
525!-- concentrations of chemical species                                   
526    IF ( air_chemistry )  THEN
527       lsp_usr = 1
528!
529!--    Chemical reactions and deposition
530       IF ( chem_gasphase_on ) THEN
531!
532!--       If required, calculate photolysis frequencies -
533!--       UNFINISHED: Why not before the intermediate timestep loop?
534          IF ( intermediate_timestep_count ==  1 )  THEN
535             CALL photolysis_control
536          ENDIF
537
538          IF ( intermediate_timestep_count == 1 .OR.                        &
539             call_chem_at_all_substeps )  THEN
540
541             CALL cpu_log( log_point_s(19), 'chem.reactions', 'start' )
542             !$OMP PARALLEL PRIVATE (i,j)
543             !$OMP DO schedule(static,1)
544             DO  i = nxl, nxr
545                DO  j = nys, nyn
546                   CALL chem_integrate (i,j)
547                ENDDO
548             ENDDO
549             !$OMP END PARALLEL
550             CALL cpu_log( log_point_s(19), 'chem.reactions', 'stop' )
551
552             IF ( deposition_dry )  THEN
553                CALL cpu_log( log_point_s(24), 'chem.deposition', 'start' )
554                DO  i = nxl, nxr
555                   DO  j = nys, nyn
556                      CALL chem_depo(i,j)
557                   ENDDO
558                ENDDO
559                CALL cpu_log( log_point_s(24), 'chem.deposition', 'stop' )
560             ENDIF
561          ENDIF
562       ENDIF
563!
564!--    Loop over chemical species       
565       CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'start' )
566       DO  lsp = 1, nspec
567          CALL exchange_horiz( chem_species(lsp)%conc, nbgp )   
568          lsp_usr = 1 
569          DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
570             IF ( TRIM(chem_species(lsp)%name) == TRIM(cs_name(lsp_usr)) )  THEN
571
572                CALL chem_boundary_conds( chem_species(lsp)%conc_p,                                 &
573                                          chem_species(lsp)%conc_pr_init ) 
574             
575             ENDIF
576             lsp_usr = lsp_usr +1
577          ENDDO
578
579
580       ENDDO
581       CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'stop' )
582
583    ENDIF       
584!
585!-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time
586!-- step. The exchange of ghost points is required after this update of the
587!-- concentrations of aerosol number and mass
588    IF ( salsa )  THEN
589       
590       IF ( time_since_reference_point >= skip_time_do_salsa )  THEN                 
591          IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  &
592          THEN
593             CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
594             !$OMP PARALLEL PRIVATE (i,j,b,c,g)
595             !$OMP DO
596!
597!--          Call salsa processes
598             DO  i = nxl, nxr
599                DO  j = nys, nyn
600                   CALL salsa_diagnostics( i, j )
601                   CALL salsa_driver( i, j, 3 )
602                   CALL salsa_diagnostics( i, j )
603                ENDDO
604             ENDDO
605             
606             CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
607             
608             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
609!
610!--          Exchange ghost points and decycle if needed.
611             DO  b = 1, nbins
612                CALL exchange_horiz( aerosol_number(b)%conc, nbgp )
613                CALL salsa_boundary_conds( aerosol_number(b)%conc_p,           &
614                                           aerosol_number(b)%init )
615                DO  c = 1, ncc_tot
616                   CALL exchange_horiz( aerosol_mass((c-1)*nbins+b)%conc, nbgp )
617                   CALL salsa_boundary_conds(                                  &
618                                           aerosol_mass((c-1)*nbins+b)%conc_p, &
619                                           aerosol_mass((c-1)*nbins+b)%init )
620                ENDDO
621             ENDDO
622             
623             IF ( .NOT. salsa_gases_from_chem )  THEN
624                DO  g = 1, ngast
625                   CALL exchange_horiz( salsa_gas(g)%conc, nbgp )
626                   CALL salsa_boundary_conds( salsa_gas(g)%conc_p,             &
627                                              salsa_gas(g)%init )
628                ENDDO
629             ENDIF
630             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
631             
632             !$OMP END PARALLEL
633             last_salsa_time = time_since_reference_point
634             
635          ENDIF
636         
637       ENDIF
638       
639    ENDIF 
640
641!
642!-- If required, calculate cloud microphysics
643    IF ( bulk_cloud_model  .AND.  .NOT. microphysics_sat_adjust  .AND.         &
644         ( intermediate_timestep_count == 1  .OR.                              &
645           call_microphysics_at_all_substeps ) )                               &
646    THEN
647       !$OMP PARALLEL PRIVATE (i,j)
648       !$OMP DO
649       DO  i = nxlg, nxrg
650          DO  j = nysg, nyng
651             CALL bcm_actions( i, j )
652           ENDDO
653       ENDDO
654       !$OMP END PARALLEL
655    ENDIF
656
657!
658!-- Loop over all prognostic equations
659!-- b, c ang g added for SALSA
660    !$OMP PARALLEL PRIVATE (i,i_omp_start,j,k,loop_start,tn,b,c,g)
661
662    !$  tn = omp_get_thread_num()
663    loop_start = .TRUE.
664
665    !$OMP DO
666    DO  i = nxl, nxr
667
668!
669!--    Store the first loop index. It differs for each thread and is required
670!--    later in advec_ws
671       IF ( loop_start )  THEN
672          loop_start  = .FALSE.
673          i_omp_start = i
674       ENDIF
675
676       DO  j = nys, nyn
677!
678!--       Tendency terms for u-velocity component. Please note, in case of
679!--       non-cyclic boundary conditions the grid point i=0 is excluded from
680!--       the prognostic equations for the u-component.   
681          IF ( i >= nxlu )  THEN
682
683             tend(:,j,i) = 0.0_wp
684             IF ( timestep_scheme(1:5) == 'runge' )  THEN
685                IF ( ws_scheme_mom )  THEN
686                   CALL advec_u_ws( i, j, i_omp_start, tn )
687                ELSE
688                   CALL advec_u_pw( i, j )
689                ENDIF
690             ELSE
691                CALL advec_u_up( i, j )
692             ENDIF
693             CALL diffusion_u( i, j )
694             CALL coriolis( i, j, 1 )
695             IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
696                CALL buoyancy( i, j, pt, 1 )
697             ENDIF
698
699!
700!--          Drag by plant canopy
701             IF ( plant_canopy )  CALL pcm_tendency( i, j, 1 )
702
703!
704!--          External pressure gradient
705             IF ( dp_external )  THEN
706                DO  k = dp_level_ind_b+1, nzt
707                   tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
708                ENDDO
709             ENDIF
710
711!
712!--          Nudging
713             IF ( nudging )  CALL nudge( i, j, simulated_time, 'u' )
714
715!
716!--          Effect of Stokes drift (in ocean mode only)
717             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 1 )
718
719!
720!--          Forces by wind turbines
721             IF ( wind_turbine )  CALL wtm_tendencies( i, j, 1 )
722
723             CALL module_interface_actions( i, j, 'u-tendency' )
724!
725!--          Prognostic equation for u-velocity component
726             DO  k = nzb+1, nzt
727
728                u_p(k,j,i) = u(k,j,i) + ( dt_3d *                               &
729                                            ( tsc(2) * tend(k,j,i) +            &
730                                              tsc(3) * tu_m(k,j,i) )            &
731                                            - tsc(5) * rdf(k)                   &
732                                                     * ( u(k,j,i) - u_init(k) ) &
733                                        ) * MERGE( 1.0_wp, 0.0_wp,              &
734                                                   BTEST( wall_flags_0(k,j,i), 1 )&
735                                                 )
736             ENDDO
737
738!
739!--          Add turbulence generated by wave breaking (in ocean mode only)
740             IF ( wave_breaking  .AND.                                         &
741               intermediate_timestep_count == intermediate_timestep_count_max )&
742             THEN
743                CALL wave_breaking_term( i, j, 1 )
744             ENDIF
745
746!
747!--          Calculate tendencies for the next Runge-Kutta step
748             IF ( timestep_scheme(1:5) == 'runge' )  THEN
749                IF ( intermediate_timestep_count == 1 )  THEN
750                   DO  k = nzb+1, nzt
751                      tu_m(k,j,i) = tend(k,j,i)
752                   ENDDO
753                ELSEIF ( intermediate_timestep_count < &
754                         intermediate_timestep_count_max )  THEN
755                   DO  k = nzb+1, nzt
756                      tu_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                &
757                                     + 5.3125_wp * tu_m(k,j,i)
758                   ENDDO
759                ENDIF
760             ENDIF
761
762          ENDIF
763!
764!--       Tendency terms for v-velocity component. Please note, in case of
765!--       non-cyclic boundary conditions the grid point j=0 is excluded from
766!--       the prognostic equations for the v-component. !--       
767          IF ( j >= nysv )  THEN
768
769             tend(:,j,i) = 0.0_wp
770             IF ( timestep_scheme(1:5) == 'runge' )  THEN
771                IF ( ws_scheme_mom )  THEN
772                    CALL advec_v_ws( i, j, i_omp_start, tn )
773                ELSE
774                    CALL advec_v_pw( i, j )
775                ENDIF
776             ELSE
777                CALL advec_v_up( i, j )
778             ENDIF
779             CALL diffusion_v( i, j )
780             CALL coriolis( i, j, 2 )
781
782!
783!--          Drag by plant canopy
784             IF ( plant_canopy )  CALL pcm_tendency( i, j, 2 )
785
786!
787!--          External pressure gradient
788             IF ( dp_external )  THEN
789                DO  k = dp_level_ind_b+1, nzt
790                   tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
791                ENDDO
792             ENDIF
793
794!
795!--          Nudging
796             IF ( nudging )  CALL nudge( i, j, simulated_time, 'v' )
797
798!
799!--          Effect of Stokes drift (in ocean mode only)
800             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 2 )
801
802!
803!--          Forces by wind turbines
804             IF ( wind_turbine )  CALL wtm_tendencies( i, j, 2 )
805
806             CALL module_interface_actions( i, j, 'v-tendency' )
807!
808!--          Prognostic equation for v-velocity component
809             DO  k = nzb+1, nzt
810                v_p(k,j,i) = v(k,j,i) + ( dt_3d *                              &
811                                            ( tsc(2) * tend(k,j,i) +           &
812                                              tsc(3) * tv_m(k,j,i) )           &
813                                            - tsc(5) * rdf(k)                  &
814                                                     * ( v(k,j,i) - v_init(k) )&
815                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
816                                                   BTEST( wall_flags_0(k,j,i), 2 )&
817                                                 )
818             ENDDO
819
820!
821!--          Add turbulence generated by wave breaking (in ocean mode only)
822             IF ( wave_breaking  .AND.                                         &
823               intermediate_timestep_count == intermediate_timestep_count_max )&
824             THEN
825                CALL wave_breaking_term( i, j, 2 )
826             ENDIF
827
828!
829!--          Calculate tendencies for the next Runge-Kutta step
830             IF ( timestep_scheme(1:5) == 'runge' )  THEN
831                IF ( intermediate_timestep_count == 1 )  THEN
832                   DO  k = nzb+1, nzt
833                      tv_m(k,j,i) = tend(k,j,i)
834                   ENDDO
835                ELSEIF ( intermediate_timestep_count < &
836                         intermediate_timestep_count_max )  THEN
837                   DO  k = nzb+1, nzt
838                      tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
839                                     + 5.3125_wp * tv_m(k,j,i)
840                   ENDDO
841                ENDIF
842             ENDIF
843
844          ENDIF
845
846!
847!--       Tendency terms for w-velocity component
848          tend(:,j,i) = 0.0_wp
849          IF ( timestep_scheme(1:5) == 'runge' )  THEN
850             IF ( ws_scheme_mom )  THEN
851                CALL advec_w_ws( i, j, i_omp_start, tn )
852             ELSE
853                CALL advec_w_pw( i, j )
854             END IF
855          ELSE
856             CALL advec_w_up( i, j )
857          ENDIF
858          CALL diffusion_w( i, j )
859          CALL coriolis( i, j, 3 )
860
861          IF ( .NOT. neutral )  THEN
862             IF ( ocean_mode )  THEN
863                CALL buoyancy( i, j, rho_ocean, 3 )
864             ELSE
865                IF ( .NOT. humidity )  THEN
866                   CALL buoyancy( i, j, pt, 3 )
867                ELSE
868                   CALL buoyancy( i, j, vpt, 3 )
869                ENDIF
870             ENDIF
871          ENDIF
872
873!
874!--       Drag by plant canopy
875          IF ( plant_canopy )  CALL pcm_tendency( i, j, 3 )
876
877!
878!--       Effect of Stokes drift (in ocean mode only)
879          IF ( stokes_force )  CALL stokes_drift_terms( i, j, 3 )
880
881!
882!--       Forces by wind turbines
883          IF ( wind_turbine )  CALL wtm_tendencies( i, j, 3 )
884
885          CALL module_interface_actions( i, j, 'w-tendency' )
886!
887!--       Prognostic equation for w-velocity component
888          DO  k = nzb+1, nzt-1
889             w_p(k,j,i) = w(k,j,i) + ( dt_3d *                                 &
890                                             ( tsc(2) * tend(k,j,i) +          &
891                                               tsc(3) * tw_m(k,j,i) )          &
892                                             - tsc(5) * rdf(k) * w(k,j,i)      &
893                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
894                                                BTEST( wall_flags_0(k,j,i), 3 )&
895                                              )
896          ENDDO
897
898!
899!--       Calculate tendencies for the next Runge-Kutta step
900          IF ( timestep_scheme(1:5) == 'runge' )  THEN
901             IF ( intermediate_timestep_count == 1 )  THEN
902                DO  k = nzb+1, nzt-1
903                   tw_m(k,j,i) = tend(k,j,i)
904                ENDDO
905             ELSEIF ( intermediate_timestep_count < &
906                      intermediate_timestep_count_max )  THEN
907                DO  k = nzb+1, nzt-1
908                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
909                                  + 5.3125_wp * tw_m(k,j,i)
910                ENDDO
911             ENDIF
912          ENDIF
913
914!
915!--       If required, compute prognostic equation for potential temperature
916          IF ( .NOT. neutral )  THEN
917!
918!--          Tendency terms for potential temperature
919             tend(:,j,i) = 0.0_wp
920             IF ( timestep_scheme(1:5) == 'runge' )  THEN
921                   IF ( ws_scheme_sca )  THEN
922                      CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, diss_s_pt,   &
923                                       flux_l_pt, diss_l_pt, i_omp_start, tn )
924                   ELSE
925                      CALL advec_s_pw( i, j, pt )
926                   ENDIF
927             ELSE
928                CALL advec_s_up( i, j, pt )
929             ENDIF
930             CALL diffusion_s( i, j, pt,                                       &
931                               surf_def_h(0)%shf, surf_def_h(1)%shf,           &
932                               surf_def_h(2)%shf,                              &
933                               surf_lsm_h%shf,    surf_usm_h%shf,              &
934                               surf_def_v(0)%shf, surf_def_v(1)%shf,           &
935                               surf_def_v(2)%shf, surf_def_v(3)%shf,           &
936                               surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,           &
937                               surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,           &
938                               surf_usm_v(0)%shf, surf_usm_v(1)%shf,           &
939                               surf_usm_v(2)%shf, surf_usm_v(3)%shf )
940
941!
942!--          Consideration of heat sources within the plant canopy
943             IF ( plant_canopy  .AND.                                          &
944                (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
945                CALL pcm_tendency( i, j, 4 )
946             ENDIF
947
948!
949!--          Large scale advection
950             IF ( large_scale_forcing )  THEN
951                CALL ls_advec( i, j, simulated_time, 'pt' )
952             ENDIF
953
954!
955!--          Nudging
956             IF ( nudging )  CALL nudge( i, j, simulated_time, 'pt' )
957
958!
959!--          If required, compute effect of large-scale subsidence/ascent
960             IF ( large_scale_subsidence  .AND.                                &
961                  .NOT. use_subsidence_tendencies )  THEN
962                CALL subsidence( i, j, tend, pt, pt_init, 2 )
963             ENDIF
964
965#if defined( __rrtmg )
966!
967!--          If required, add tendency due to radiative heating/cooling
968             IF ( radiation  .AND.                                             &
969                  simulated_time > skip_time_do_radiation )  THEN
970                CALL radiation_tendency ( i, j, tend )
971             ENDIF
972#endif
973
974             CALL module_interface_actions( i, j, 'pt-tendency' )
975!
976!--          Prognostic equation for potential temperature
977             DO  k = nzb+1, nzt
978                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d *                            &
979                                                  ( tsc(2) * tend(k,j,i) +     &
980                                                    tsc(3) * tpt_m(k,j,i) )    &
981                                                  - tsc(5)                     &
982                                                  * ( pt(k,j,i) - pt_init(k) ) &
983                                                  * ( rdf_sc(k) + ptdf_x(i)    &
984                                                                + ptdf_y(j) )  &
985                                          )                                    &
986                                       * MERGE( 1.0_wp, 0.0_wp,                &
987                                                BTEST( wall_flags_0(k,j,i), 0 )&
988                                              )
989             ENDDO
990
991!
992!--          Calculate tendencies for the next Runge-Kutta step
993             IF ( timestep_scheme(1:5) == 'runge' )  THEN
994                IF ( intermediate_timestep_count == 1 )  THEN
995                   DO  k = nzb+1, nzt
996                      tpt_m(k,j,i) = tend(k,j,i)
997                   ENDDO
998                ELSEIF ( intermediate_timestep_count < &
999                         intermediate_timestep_count_max )  THEN
1000                   DO  k = nzb+1, nzt
1001                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
1002                                        5.3125_wp * tpt_m(k,j,i)
1003                   ENDDO
1004                ENDIF
1005             ENDIF
1006
1007          ENDIF
1008
1009!
1010!--       If required, compute prognostic equation for total water content
1011          IF ( humidity )  THEN
1012
1013!
1014!--          Tendency-terms for total water content / scalar
1015             tend(:,j,i) = 0.0_wp
1016             IF ( timestep_scheme(1:5) == 'runge' ) &
1017             THEN
1018                IF ( ws_scheme_sca )  THEN
1019                   CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
1020                                diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
1021                ELSE
1022                   CALL advec_s_pw( i, j, q )
1023                ENDIF
1024             ELSE
1025                CALL advec_s_up( i, j, q )
1026             ENDIF
1027             CALL diffusion_s( i, j, q,                                        &
1028                               surf_def_h(0)%qsws, surf_def_h(1)%qsws,         &
1029                               surf_def_h(2)%qsws,                             &
1030                               surf_lsm_h%qsws,    surf_usm_h%qsws,            &
1031                               surf_def_v(0)%qsws, surf_def_v(1)%qsws,         &
1032                               surf_def_v(2)%qsws, surf_def_v(3)%qsws,         &
1033                               surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,         &
1034                               surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,         &
1035                               surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,         &
1036                               surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
1037
1038!
1039!--          Sink or source of humidity due to canopy elements
1040             IF ( plant_canopy )  CALL pcm_tendency( i, j, 5 )
1041
1042!
1043!--          Large scale advection
1044             IF ( large_scale_forcing )  THEN
1045                CALL ls_advec( i, j, simulated_time, 'q' )
1046             ENDIF
1047
1048!
1049!--          Nudging
1050             IF ( nudging )  CALL nudge( i, j, simulated_time, 'q' )
1051
1052!
1053!--          If required compute influence of large-scale subsidence/ascent
1054             IF ( large_scale_subsidence  .AND.                                &
1055                  .NOT. use_subsidence_tendencies )  THEN
1056                CALL subsidence( i, j, tend, q, q_init, 3 )
1057             ENDIF
1058
1059             CALL module_interface_actions( i, j, 'q-tendency' )
1060
1061!
1062!--          Prognostic equation for total water content / scalar
1063             DO  k = nzb+1, nzt
1064                q_p(k,j,i) = q(k,j,i) + ( dt_3d *                              &
1065                                                ( tsc(2) * tend(k,j,i) +       &
1066                                                  tsc(3) * tq_m(k,j,i) )       &
1067                                                - tsc(5) * rdf_sc(k) *         &
1068                                                      ( q(k,j,i) - q_init(k) ) &
1069                                        )                                      &
1070                                       * MERGE( 1.0_wp, 0.0_wp,                &
1071                                                BTEST( wall_flags_0(k,j,i), 0 )&
1072                                              )               
1073                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
1074             ENDDO
1075
1076!
1077!--          Calculate tendencies for the next Runge-Kutta step
1078             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1079                IF ( intermediate_timestep_count == 1 )  THEN
1080                   DO  k = nzb+1, nzt
1081                      tq_m(k,j,i) = tend(k,j,i)
1082                   ENDDO
1083                ELSEIF ( intermediate_timestep_count < &
1084                         intermediate_timestep_count_max )  THEN
1085                   DO  k = nzb+1, nzt
1086                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
1087                                       5.3125_wp * tq_m(k,j,i)
1088                   ENDDO
1089                ENDIF
1090             ENDIF
1091
1092!
1093!--          If required, calculate prognostic equations for cloud water content
1094!--          and cloud drop concentration
1095             IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
1096!
1097!--             Calculate prognostic equation for cloud water content
1098                tend(:,j,i) = 0.0_wp
1099                IF ( timestep_scheme(1:5) == 'runge' ) &
1100                THEN
1101                   IF ( ws_scheme_sca )  THEN
1102                      CALL advec_s_ws( i, j, qc, 'qc', flux_s_qc,       &
1103                                       diss_s_qc, flux_l_qc, diss_l_qc, &
1104                                       i_omp_start, tn )
1105                   ELSE
1106                      CALL advec_s_pw( i, j, qc )
1107                   ENDIF
1108                ELSE
1109                   CALL advec_s_up( i, j, qc )
1110                ENDIF
1111                CALL diffusion_s( i, j, qc,                                   &
1112                                  surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,   &
1113                                  surf_def_h(2)%qcsws,                        &
1114                                  surf_lsm_h%qcsws,    surf_usm_h%qcsws,      & 
1115                                  surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,   &
1116                                  surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,   &
1117                                  surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,   &
1118                                  surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,   &
1119                                  surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,   &
1120                                  surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
1121
1122!
1123!--             Prognostic equation for cloud water content
1124                DO  k = nzb+1, nzt
1125                   qc_p(k,j,i) = qc(k,j,i) + ( dt_3d *                         &
1126                                                      ( tsc(2) * tend(k,j,i) + &
1127                                                        tsc(3) * tqc_m(k,j,i) )&
1128                                                      - tsc(5) * rdf_sc(k)     &
1129                                                               * qc(k,j,i)     &
1130                                             )                                 &
1131                                       * MERGE( 1.0_wp, 0.0_wp,                &
1132                                                BTEST( wall_flags_0(k,j,i), 0 )&
1133                                              ) 
1134                   IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
1135                ENDDO
1136!
1137!--             Calculate tendencies for the next Runge-Kutta step
1138                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1139                   IF ( intermediate_timestep_count == 1 )  THEN
1140                      DO  k = nzb+1, nzt
1141                         tqc_m(k,j,i) = tend(k,j,i)
1142                      ENDDO
1143                   ELSEIF ( intermediate_timestep_count < &
1144                            intermediate_timestep_count_max )  THEN
1145                      DO  k = nzb+1, nzt
1146                         tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
1147                                           5.3125_wp * tqc_m(k,j,i)
1148                      ENDDO
1149                   ENDIF
1150                ENDIF
1151
1152!
1153!--             Calculate prognostic equation for cloud drop concentration.
1154                tend(:,j,i) = 0.0_wp
1155                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1156                   IF ( ws_scheme_sca )  THEN
1157                      CALL advec_s_ws( i, j, nc, 'nc', flux_s_nc,    &
1158                                    diss_s_nc, flux_l_nc, diss_l_nc, &
1159                                    i_omp_start, tn )
1160                   ELSE
1161                      CALL advec_s_pw( i, j, nc )
1162                   ENDIF
1163                ELSE
1164                   CALL advec_s_up( i, j, nc )
1165                ENDIF
1166                CALL diffusion_s( i, j, nc,                                    &
1167                                  surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,    &
1168                                  surf_def_h(2)%ncsws,                         &
1169                                  surf_lsm_h%ncsws,    surf_usm_h%ncsws,       &
1170                                  surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,    &
1171                                  surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,    &
1172                                  surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,    &
1173                                  surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,    &
1174                                  surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,    &
1175                                  surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
1176
1177!
1178!--             Prognostic equation for cloud drop concentration
1179                DO  k = nzb+1, nzt
1180                   nc_p(k,j,i) = nc(k,j,i) + ( dt_3d *                         &
1181                                                      ( tsc(2) * tend(k,j,i) + &
1182                                                        tsc(3) * tnc_m(k,j,i) )&
1183                                                      - tsc(5) * rdf_sc(k)     &
1184                                                               * nc(k,j,i)     &
1185                                             )                                 &
1186                                       * MERGE( 1.0_wp, 0.0_wp,                &
1187                                                BTEST( wall_flags_0(k,j,i), 0 )&
1188                                              )
1189                   IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
1190                ENDDO
1191!
1192!--             Calculate tendencies for the next Runge-Kutta step
1193                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1194                   IF ( intermediate_timestep_count == 1 )  THEN
1195                      DO  k = nzb+1, nzt
1196                         tnc_m(k,j,i) = tend(k,j,i)
1197                      ENDDO
1198                   ELSEIF ( intermediate_timestep_count < &
1199                            intermediate_timestep_count_max )  THEN
1200                      DO  k = nzb+1, nzt
1201                         tnc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
1202                                           5.3125_wp * tnc_m(k,j,i)
1203                      ENDDO
1204                   ENDIF
1205                ENDIF
1206
1207             ENDIF
1208!
1209!--          If required, calculate prognostic equations for rain water content
1210!--          and rain drop concentration
1211             IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
1212!
1213!--             Calculate prognostic equation for rain water content
1214                tend(:,j,i) = 0.0_wp
1215                IF ( timestep_scheme(1:5) == 'runge' ) &
1216                THEN
1217                   IF ( ws_scheme_sca )  THEN
1218                      CALL advec_s_ws( i, j, qr, 'qr', flux_s_qr,       &
1219                                       diss_s_qr, flux_l_qr, diss_l_qr, &
1220                                       i_omp_start, tn )
1221                   ELSE
1222                      CALL advec_s_pw( i, j, qr )
1223                   ENDIF
1224                ELSE
1225                   CALL advec_s_up( i, j, qr )
1226                ENDIF
1227                CALL diffusion_s( i, j, qr,                                   &
1228                                  surf_def_h(0)%qrsws, surf_def_h(1)%qrsws,   &
1229                                  surf_def_h(2)%qrsws,                        &
1230                                  surf_lsm_h%qrsws,    surf_usm_h%qrsws,      & 
1231                                  surf_def_v(0)%qrsws, surf_def_v(1)%qrsws,   &
1232                                  surf_def_v(2)%qrsws, surf_def_v(3)%qrsws,   &
1233                                  surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws,   &
1234                                  surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws,   &
1235                                  surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws,   &
1236                                  surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws )
1237
1238!
1239!--             Prognostic equation for rain water content
1240                DO  k = nzb+1, nzt
1241                   qr_p(k,j,i) = qr(k,j,i) + ( dt_3d *                         &
1242                                                      ( tsc(2) * tend(k,j,i) + &
1243                                                        tsc(3) * tqr_m(k,j,i) )&
1244                                                      - tsc(5) * rdf_sc(k)     &
1245                                                               * qr(k,j,i)     &
1246                                             )                                 &
1247                                       * MERGE( 1.0_wp, 0.0_wp,                &
1248                                                BTEST( wall_flags_0(k,j,i), 0 )&
1249                                              ) 
1250                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
1251                ENDDO
1252!
1253!--             Calculate tendencies for the next Runge-Kutta step
1254                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1255                   IF ( intermediate_timestep_count == 1 )  THEN
1256                      DO  k = nzb+1, nzt
1257                         tqr_m(k,j,i) = tend(k,j,i)
1258                      ENDDO
1259                   ELSEIF ( intermediate_timestep_count < &
1260                            intermediate_timestep_count_max )  THEN
1261                      DO  k = nzb+1, nzt
1262                         tqr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
1263                                           5.3125_wp * tqr_m(k,j,i)
1264                      ENDDO
1265                   ENDIF
1266                ENDIF
1267
1268!
1269!--             Calculate prognostic equation for rain drop concentration.
1270                tend(:,j,i) = 0.0_wp
1271                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1272                   IF ( ws_scheme_sca )  THEN
1273                      CALL advec_s_ws( i, j, nr, 'nr', flux_s_nr,    &
1274                                    diss_s_nr, flux_l_nr, diss_l_nr, &
1275                                    i_omp_start, tn )
1276                   ELSE
1277                      CALL advec_s_pw( i, j, nr )
1278                   ENDIF
1279                ELSE
1280                   CALL advec_s_up( i, j, nr )
1281                ENDIF
1282                CALL diffusion_s( i, j, nr,                                    &
1283                                  surf_def_h(0)%nrsws, surf_def_h(1)%nrsws,    &
1284                                  surf_def_h(2)%nrsws,                         &
1285                                  surf_lsm_h%nrsws,    surf_usm_h%nrsws,       &
1286                                  surf_def_v(0)%nrsws, surf_def_v(1)%nrsws,    &
1287                                  surf_def_v(2)%nrsws, surf_def_v(3)%nrsws,    &
1288                                  surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws,    &
1289                                  surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws,    &
1290                                  surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws,    &
1291                                  surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws )
1292
1293!
1294!--             Prognostic equation for rain drop concentration
1295                DO  k = nzb+1, nzt
1296                   nr_p(k,j,i) = nr(k,j,i) + ( dt_3d *                         &
1297                                                      ( tsc(2) * tend(k,j,i) + &
1298                                                        tsc(3) * tnr_m(k,j,i) )&
1299                                                      - tsc(5) * rdf_sc(k)     &
1300                                                               * nr(k,j,i)     &
1301                                             )                                 &
1302                                       * MERGE( 1.0_wp, 0.0_wp,                &
1303                                                BTEST( wall_flags_0(k,j,i), 0 )&
1304                                              )
1305                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
1306                ENDDO
1307!
1308!--             Calculate tendencies for the next Runge-Kutta step
1309                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1310                   IF ( intermediate_timestep_count == 1 )  THEN
1311                      DO  k = nzb+1, nzt
1312                         tnr_m(k,j,i) = tend(k,j,i)
1313                      ENDDO
1314                   ELSEIF ( intermediate_timestep_count < &
1315                            intermediate_timestep_count_max )  THEN
1316                      DO  k = nzb+1, nzt
1317                         tnr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
1318                                           5.3125_wp * tnr_m(k,j,i)
1319                      ENDDO
1320                   ENDIF
1321                ENDIF
1322
1323             ENDIF
1324
1325          ENDIF
1326
1327!
1328!--       If required, compute prognostic equation for scalar
1329          IF ( passive_scalar )  THEN
1330!
1331!--          Tendency-terms for total water content / scalar
1332             tend(:,j,i) = 0.0_wp
1333             IF ( timestep_scheme(1:5) == 'runge' ) &
1334             THEN
1335                IF ( ws_scheme_sca )  THEN
1336                   CALL advec_s_ws( i, j, s, 's', flux_s_s, &
1337                                diss_s_s, flux_l_s, diss_l_s, i_omp_start, tn )
1338                ELSE
1339                   CALL advec_s_pw( i, j, s )
1340                ENDIF
1341             ELSE
1342                CALL advec_s_up( i, j, s )
1343             ENDIF
1344             CALL diffusion_s( i, j, s,                                        &
1345                               surf_def_h(0)%ssws, surf_def_h(1)%ssws,         &
1346                               surf_def_h(2)%ssws,                             &
1347                               surf_lsm_h%ssws,    surf_usm_h%ssws,            &
1348                               surf_def_v(0)%ssws, surf_def_v(1)%ssws,         &
1349                               surf_def_v(2)%ssws, surf_def_v(3)%ssws,         &
1350                               surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,         &
1351                               surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,         &
1352                               surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,         &
1353                               surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1354
1355!
1356!--          Sink or source of scalar concentration due to canopy elements
1357             IF ( plant_canopy )  CALL pcm_tendency( i, j, 7 )
1358
1359!
1360!--          Large scale advection, still need to be extended for scalars
1361!              IF ( large_scale_forcing )  THEN
1362!                 CALL ls_advec( i, j, simulated_time, 's' )
1363!              ENDIF
1364
1365!
1366!--          Nudging, still need to be extended for scalars
1367!              IF ( nudging )  CALL nudge( i, j, simulated_time, 's' )
1368
1369!
1370!--          If required compute influence of large-scale subsidence/ascent.
1371!--          Note, the last argument is of no meaning in this case, as it is
1372!--          only used in conjunction with large_scale_forcing, which is to
1373!--          date not implemented for scalars.
1374             IF ( large_scale_subsidence  .AND.                                &
1375                  .NOT. use_subsidence_tendencies  .AND.                       &
1376                  .NOT. large_scale_forcing )  THEN
1377                CALL subsidence( i, j, tend, s, s_init, 3 )
1378             ENDIF
1379
1380             CALL module_interface_actions( i, j, 's-tendency' )
1381
1382!
1383!--          Prognostic equation for scalar
1384             DO  k = nzb+1, nzt
1385                s_p(k,j,i) = s(k,j,i) + (  dt_3d *                             &
1386                                            ( tsc(2) * tend(k,j,i) +           &
1387                                              tsc(3) * ts_m(k,j,i) )           &
1388                                            - tsc(5) * rdf_sc(k)               &
1389                                                     * ( s(k,j,i) - s_init(k) )&
1390                                        )                                      &
1391                                       * MERGE( 1.0_wp, 0.0_wp,                &
1392                                                BTEST( wall_flags_0(k,j,i), 0 )&
1393                                              )
1394                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
1395             ENDDO
1396
1397!
1398!--          Calculate tendencies for the next Runge-Kutta step
1399             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1400                IF ( intermediate_timestep_count == 1 )  THEN
1401                   DO  k = nzb+1, nzt
1402                      ts_m(k,j,i) = tend(k,j,i)
1403                   ENDDO
1404                ELSEIF ( intermediate_timestep_count < &
1405                         intermediate_timestep_count_max )  THEN
1406                   DO  k = nzb+1, nzt
1407                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
1408                                       5.3125_wp * ts_m(k,j,i)
1409                   ENDDO
1410                ENDIF
1411             ENDIF
1412
1413          ENDIF
1414!
1415!--       Calculate prognostic equations for turbulence closure
1416          CALL tcm_prognostic_equations( i, j, i_omp_start, tn )
1417
1418!
1419!--       Calculate prognostic equation for chemical quantites
1420          IF ( air_chemistry )  THEN
1421             !> TODO: remove time measurement since it slows down performance because it will be called extremely often
1422!
1423!--          Loop over chemical species
1424             DO  lsp = 1, nvar                         
1425                CALL chem_prognostic_equations( chem_species(lsp)%conc_p,      &
1426                                     chem_species(lsp)%conc,                   &
1427                                     chem_species(lsp)%tconc_m,                &
1428                                     chem_species(lsp)%conc_pr_init,           &
1429                                     i, j, i_omp_start, tn, lsp,               &
1430                                     chem_species(lsp)%flux_s_cs,              &
1431                                     chem_species(lsp)%diss_s_cs,              &
1432                                     chem_species(lsp)%flux_l_cs,              &
1433                                     chem_species(lsp)%diss_l_cs )       
1434             ENDDO
1435         
1436          ENDIF   ! Chemical equations
1437!
1438!--       Calculate prognostic equations for the ocean
1439          IF ( ocean_mode )  THEN
1440             CALL ocean_prognostic_equations( i, j, i_omp_start, tn )
1441          ENDIF
1442       
1443          IF ( salsa )  THEN 
1444!
1445!--          Loop over aerosol size bins: number and mass bins
1446             IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
1447             
1448                DO  b = 1, nbins               
1449                   sums_salsa_ws_l = aerosol_number(b)%sums_ws_l
1450                   CALL salsa_tendency( 'aerosol_number',                      &
1451                                         aerosol_number(b)%conc_p,             &
1452                                         aerosol_number(b)%conc,               &
1453                                         aerosol_number(b)%tconc_m,            &
1454                                         i, j, i_omp_start, tn, b, b,          &
1455                                         aerosol_number(b)%flux_s,             &
1456                                         aerosol_number(b)%diss_s,             &
1457                                         aerosol_number(b)%flux_l,             &
1458                                         aerosol_number(b)%diss_l,             &
1459                                         aerosol_number(b)%init )
1460                   aerosol_number(b)%sums_ws_l = sums_salsa_ws_l
1461                   DO  c = 1, ncc_tot
1462                      sums_salsa_ws_l = aerosol_mass((c-1)*nbins+b)%sums_ws_l 
1463                      CALL salsa_tendency( 'aerosol_mass',                     &
1464                                            aerosol_mass((c-1)*nbins+b)%conc_p,&
1465                                            aerosol_mass((c-1)*nbins+b)%conc,  &
1466                                            aerosol_mass((c-1)*nbins+b)%tconc_m,&
1467                                            i, j, i_omp_start, tn, b, c,       &
1468                                            aerosol_mass((c-1)*nbins+b)%flux_s,&
1469                                            aerosol_mass((c-1)*nbins+b)%diss_s,&
1470                                            aerosol_mass((c-1)*nbins+b)%flux_l,&
1471                                            aerosol_mass((c-1)*nbins+b)%diss_l,&
1472                                            aerosol_mass((c-1)*nbins+b)%init )
1473                      aerosol_mass((c-1)*nbins+b)%sums_ws_l = sums_salsa_ws_l
1474                   ENDDO
1475                ENDDO
1476                IF ( .NOT. salsa_gases_from_chem )  THEN
1477                   DO  g = 1, ngast
1478                      sums_salsa_ws_l = salsa_gas(g)%sums_ws_l
1479                      CALL salsa_tendency( 'salsa_gas', salsa_gas(g)%conc_p,   &
1480                                      salsa_gas(g)%conc, salsa_gas(g)%tconc_m, &
1481                                      i, j, i_omp_start, tn, g, g,             &
1482                                      salsa_gas(g)%flux_s, salsa_gas(g)%diss_s,&
1483                                      salsa_gas(g)%flux_l, salsa_gas(g)%diss_l,&
1484                                      salsa_gas(g)%init )
1485                      salsa_gas(g)%sums_ws_l = sums_salsa_ws_l
1486                   ENDDO
1487                ENDIF
1488             
1489             ENDIF
1490             
1491          ENDIF       
1492
1493       ENDDO  ! loop over j
1494    ENDDO  ! loop over i
1495    !$OMP END PARALLEL
1496
1497
1498    CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1499
1500
1501 END SUBROUTINE prognostic_equations_cache
1502
1503
1504!------------------------------------------------------------------------------!
1505! Description:
1506! ------------
1507!> Version for vector machines
1508!------------------------------------------------------------------------------!
1509
1510 SUBROUTINE prognostic_equations_vector
1511
1512
1513    IMPLICIT NONE
1514
1515    INTEGER(iwp) ::  b     !< index for aerosol size bins (salsa)
1516    INTEGER(iwp) ::  c     !< index for chemical compounds (salsa)
1517    INTEGER(iwp) ::  g     !< index for gaseous compounds (salsa)
1518    INTEGER(iwp) ::  i     !<
1519    INTEGER(iwp) ::  j     !<
1520    INTEGER(iwp) ::  k     !<
1521    INTEGER(iwp) ::  lsp   !< running index for chemical species
1522
1523    REAL(wp)     ::  sbt  !<
1524
1525!
1526!-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time
1527!-- step. The exchange of ghost points is required after this update of the
1528!-- concentrations of aerosol number and mass
1529    IF ( salsa )  THEN
1530       
1531       IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
1532       
1533          IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  &
1534          THEN
1535         
1536             CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
1537             !$OMP PARALLEL PRIVATE (i,j,b,c,g)
1538             !$OMP DO
1539!
1540!--          Call salsa processes
1541             DO  i = nxl, nxr
1542                DO  j = nys, nyn
1543                   CALL salsa_diagnostics( i, j )
1544                   CALL salsa_driver( i, j, 3 )
1545                   CALL salsa_diagnostics( i, j )
1546                ENDDO
1547             ENDDO
1548             
1549             CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
1550             
1551             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
1552!
1553!--          Exchange ghost points and decycle if needed.
1554             DO  b = 1, nbins
1555                CALL exchange_horiz( aerosol_number(b)%conc, nbgp )
1556                CALL salsa_boundary_conds( aerosol_number(b)%conc_p,           &
1557                                           aerosol_number(b)%init )
1558                DO  c = 1, ncc_tot
1559                   CALL exchange_horiz( aerosol_mass((c-1)*nbins+b)%conc, nbgp )
1560                   CALL salsa_boundary_conds(                                  &
1561                                           aerosol_mass((c-1)*nbins+b)%conc_p, &
1562                                           aerosol_mass((c-1)*nbins+b)%init )
1563                ENDDO
1564             ENDDO
1565             
1566             IF ( .NOT. salsa_gases_from_chem )  THEN
1567                DO  g = 1, ngast
1568                   CALL exchange_horiz( salsa_gas(g)%conc, nbgp )
1569                   CALL salsa_boundary_conds( salsa_gas(g)%conc_p,             &
1570                                              salsa_gas(g)%init )
1571                ENDDO
1572             ENDIF
1573             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
1574             
1575             !$OMP END PARALLEL
1576             last_salsa_time = time_since_reference_point
1577             
1578          ENDIF
1579         
1580       ENDIF
1581       
1582    ENDIF 
1583
1584!
1585!-- If required, calculate cloud microphysical impacts
1586    IF ( bulk_cloud_model  .AND.  .NOT. microphysics_sat_adjust  .AND.         &
1587         ( intermediate_timestep_count == 1  .OR.                              &
1588           call_microphysics_at_all_substeps )                                 &
1589       )  THEN
1590       CALL cpu_log( log_point(51), 'microphysics', 'start' )
1591       CALL bcm_actions
1592       CALL cpu_log( log_point(51), 'microphysics', 'stop' )
1593    ENDIF
1594
1595!
1596!-- u-velocity component
1597    CALL cpu_log( log_point(5), 'u-equation', 'start' )
1598
1599    !$ACC KERNELS PRESENT(tend)
1600    tend = 0.0_wp
1601    !$ACC END KERNELS
1602    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1603       IF ( ws_scheme_mom )  THEN
1604          CALL advec_u_ws
1605       ELSE
1606          CALL advec_u_pw
1607       ENDIF
1608    ELSE
1609       CALL advec_u_up
1610    ENDIF
1611    CALL diffusion_u
1612    CALL coriolis( 1 )
1613    IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
1614       CALL buoyancy( pt, 1 )
1615    ENDIF
1616
1617!
1618!-- Drag by plant canopy
1619    IF ( plant_canopy )  CALL pcm_tendency( 1 )
1620
1621!
1622!-- External pressure gradient
1623    IF ( dp_external )  THEN
1624       DO  i = nxlu, nxr
1625          DO  j = nys, nyn
1626             DO  k = dp_level_ind_b+1, nzt
1627                tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
1628             ENDDO
1629          ENDDO
1630       ENDDO
1631    ENDIF
1632
1633!
1634!-- Nudging
1635    IF ( nudging )  CALL nudge( simulated_time, 'u' )
1636
1637!
1638!-- Effect of Stokes drift (in ocean mode only)
1639    IF ( stokes_force )  CALL stokes_drift_terms( 1 )
1640
1641!
1642!-- Forces by wind turbines
1643    IF ( wind_turbine )  CALL wtm_tendencies( 1 )
1644
1645    CALL module_interface_actions( 'u-tendency' )
1646
1647!
1648!-- Prognostic equation for u-velocity component
1649    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1650    !$ACC PRESENT(u, tend, tu_m, u_init, rdf, wall_flags_0) &
1651    !$ACC PRESENT(tsc(2:5)) &
1652    !$ACC PRESENT(u_p)
1653    DO  i = nxlu, nxr
1654       DO  j = nys, nyn
1655          DO  k = nzb+1, nzt
1656             u_p(k,j,i) = u(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +          &
1657                                                 tsc(3) * tu_m(k,j,i) )          &
1658                                               - tsc(5) * rdf(k) *               &
1659                                                        ( u(k,j,i) - u_init(k) ) &
1660                                     ) * MERGE( 1.0_wp, 0.0_wp,                  &
1661                                                 BTEST( wall_flags_0(k,j,i), 1 ) &
1662                                              )
1663          ENDDO
1664       ENDDO
1665    ENDDO
1666
1667!
1668!-- Add turbulence generated by wave breaking (in ocean mode only)
1669    IF ( wave_breaking  .AND.                                                  &
1670         intermediate_timestep_count == intermediate_timestep_count_max )      &
1671    THEN
1672       CALL wave_breaking_term( 1 )
1673    ENDIF
1674
1675!
1676!-- Calculate tendencies for the next Runge-Kutta step
1677    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1678       IF ( intermediate_timestep_count == 1 )  THEN
1679          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1680          !$ACC PRESENT(tend, tu_m)
1681          DO  i = nxlu, nxr
1682             DO  j = nys, nyn
1683                DO  k = nzb+1, nzt
1684                   tu_m(k,j,i) = tend(k,j,i)
1685                ENDDO
1686             ENDDO
1687          ENDDO
1688       ELSEIF ( intermediate_timestep_count < &
1689                intermediate_timestep_count_max )  THEN
1690          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1691          !$ACC PRESENT(tend, tu_m)
1692          DO  i = nxlu, nxr
1693             DO  j = nys, nyn
1694                DO  k = nzb+1, nzt
1695                   tu_m(k,j,i) =    -9.5625_wp * tend(k,j,i)                   &
1696                                   + 5.3125_wp * tu_m(k,j,i)
1697                ENDDO
1698             ENDDO
1699          ENDDO
1700       ENDIF
1701    ENDIF
1702
1703    CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1704
1705!
1706!-- v-velocity component
1707    CALL cpu_log( log_point(6), 'v-equation', 'start' )
1708
1709    !$ACC KERNELS PRESENT(tend)
1710    tend = 0.0_wp
1711    !$ACC END KERNELS
1712    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1713       IF ( ws_scheme_mom )  THEN
1714          CALL advec_v_ws
1715       ELSE
1716          CALL advec_v_pw
1717       END IF
1718    ELSE
1719       CALL advec_v_up
1720    ENDIF
1721    CALL diffusion_v
1722    CALL coriolis( 2 )
1723
1724!
1725!-- Drag by plant canopy
1726    IF ( plant_canopy )  CALL pcm_tendency( 2 )
1727
1728!
1729!-- External pressure gradient
1730    IF ( dp_external )  THEN
1731       DO  i = nxl, nxr
1732          DO  j = nysv, nyn
1733             DO  k = dp_level_ind_b+1, nzt
1734                tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1735             ENDDO
1736          ENDDO
1737       ENDDO
1738    ENDIF
1739
1740!
1741!-- Nudging
1742    IF ( nudging )  CALL nudge( simulated_time, 'v' )
1743
1744!
1745!-- Effect of Stokes drift (in ocean mode only)
1746    IF ( stokes_force )  CALL stokes_drift_terms( 2 )
1747
1748!
1749!-- Forces by wind turbines
1750    IF ( wind_turbine )  CALL wtm_tendencies( 2 )
1751
1752    CALL module_interface_actions( 'v-tendency' )
1753
1754!
1755!-- Prognostic equation for v-velocity component
1756    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1757    !$ACC PRESENT(v, tend, tv_m, v_init, rdf, wall_flags_0) &
1758    !$ACC PRESENT(tsc(2:5)) &
1759    !$ACC PRESENT(v_p)
1760    DO  i = nxl, nxr
1761       DO  j = nysv, nyn
1762          DO  k = nzb+1, nzt
1763             v_p(k,j,i) = v(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1764                                                 tsc(3) * tv_m(k,j,i) )        &
1765                                               - tsc(5) * rdf(k) *             &
1766                                                      ( v(k,j,i) - v_init(k) ) &
1767                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1768                                                BTEST( wall_flags_0(k,j,i), 2 )&
1769                                              )
1770          ENDDO
1771       ENDDO
1772    ENDDO
1773
1774!
1775!-- Add turbulence generated by wave breaking (in ocean mode only)
1776    IF ( wave_breaking  .AND.                                                  &
1777         intermediate_timestep_count == intermediate_timestep_count_max )      &
1778    THEN
1779       CALL wave_breaking_term( 2 )
1780    ENDIF
1781
1782!
1783!-- Calculate tendencies for the next Runge-Kutta step
1784    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1785       IF ( intermediate_timestep_count == 1 )  THEN
1786          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1787          !$ACC PRESENT(tend, tv_m)
1788          DO  i = nxl, nxr
1789             DO  j = nysv, nyn
1790                DO  k = nzb+1, nzt
1791                   tv_m(k,j,i) = tend(k,j,i)
1792                ENDDO
1793             ENDDO
1794          ENDDO
1795       ELSEIF ( intermediate_timestep_count < &
1796                intermediate_timestep_count_max )  THEN
1797          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1798          !$ACC PRESENT(tend, tv_m)
1799          DO  i = nxl, nxr
1800             DO  j = nysv, nyn
1801                DO  k = nzb+1, nzt
1802                   tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1803                                  + 5.3125_wp * tv_m(k,j,i)
1804                ENDDO
1805             ENDDO
1806          ENDDO
1807       ENDIF
1808    ENDIF
1809
1810    CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1811
1812!
1813!-- w-velocity component
1814    CALL cpu_log( log_point(7), 'w-equation', 'start' )
1815
1816    !$ACC KERNELS PRESENT(tend)
1817    tend = 0.0_wp
1818    !$ACC END KERNELS
1819    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1820       IF ( ws_scheme_mom )  THEN
1821          CALL advec_w_ws
1822       ELSE
1823          CALL advec_w_pw
1824       ENDIF
1825    ELSE
1826       CALL advec_w_up
1827    ENDIF
1828    CALL diffusion_w
1829    CALL coriolis( 3 )
1830
1831    IF ( .NOT. neutral )  THEN
1832       IF ( ocean_mode )  THEN
1833          CALL buoyancy( rho_ocean, 3 )
1834       ELSE
1835          IF ( .NOT. humidity )  THEN
1836             CALL buoyancy( pt, 3 )
1837          ELSE
1838             CALL buoyancy( vpt, 3 )
1839          ENDIF
1840       ENDIF
1841    ENDIF
1842
1843!
1844!-- Drag by plant canopy
1845    IF ( plant_canopy )  CALL pcm_tendency( 3 )
1846
1847!
1848!-- Effect of Stokes drift (in ocean mode only)
1849    IF ( stokes_force )  CALL stokes_drift_terms( 3 )
1850
1851!
1852!-- Forces by wind turbines
1853    IF ( wind_turbine )  CALL wtm_tendencies( 3 )
1854
1855    CALL module_interface_actions( 'w-tendency' )
1856
1857!
1858!-- Prognostic equation for w-velocity component
1859    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1860    !$ACC PRESENT(w, tend, tw_m, v_init, rdf, wall_flags_0) &
1861    !$ACC PRESENT(tsc(2:5)) &
1862    !$ACC PRESENT(w_p)
1863    DO  i = nxl, nxr
1864       DO  j = nys, nyn
1865          DO  k = nzb+1, nzt-1
1866             w_p(k,j,i) = w(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1867                                                 tsc(3) * tw_m(k,j,i) )        &
1868                                               - tsc(5) * rdf(k) * w(k,j,i)    &
1869                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1870                                                BTEST( wall_flags_0(k,j,i), 3 )&
1871                                              )
1872          ENDDO
1873       ENDDO
1874    ENDDO
1875
1876!
1877!-- Calculate tendencies for the next Runge-Kutta step
1878    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1879       IF ( intermediate_timestep_count == 1 )  THEN
1880          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1881          !$ACC PRESENT(tend, tw_m)
1882          DO  i = nxl, nxr
1883             DO  j = nys, nyn
1884                DO  k = nzb+1, nzt-1
1885                   tw_m(k,j,i) = tend(k,j,i)
1886                ENDDO
1887             ENDDO
1888          ENDDO
1889       ELSEIF ( intermediate_timestep_count < &
1890                intermediate_timestep_count_max )  THEN
1891          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1892          !$ACC PRESENT(tend, tw_m)
1893          DO  i = nxl, nxr
1894             DO  j = nys, nyn
1895                DO  k = nzb+1, nzt-1
1896                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1897                                  + 5.3125_wp * tw_m(k,j,i)
1898                ENDDO
1899             ENDDO
1900          ENDDO
1901       ENDIF
1902    ENDIF
1903
1904    CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1905
1906
1907!
1908!-- If required, compute prognostic equation for potential temperature
1909    IF ( .NOT. neutral )  THEN
1910
1911       CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1912
1913!
1914!--    pt-tendency terms with communication
1915       sbt = tsc(2)
1916       IF ( scalar_advec == 'bc-scheme' )  THEN
1917
1918          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1919!
1920!--          Bott-Chlond scheme always uses Euler time step. Thus:
1921             sbt = 1.0_wp
1922          ENDIF
1923          tend = 0.0_wp
1924          CALL advec_s_bc( pt, 'pt' )
1925
1926       ENDIF
1927
1928!
1929!--    pt-tendency terms with no communication
1930       IF ( scalar_advec /= 'bc-scheme' )  THEN
1931          !$ACC KERNELS PRESENT(tend)
1932          tend = 0.0_wp
1933          !$ACC END KERNELS
1934          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1935             IF ( ws_scheme_sca )  THEN
1936                CALL advec_s_ws( pt, 'pt' )
1937             ELSE
1938                CALL advec_s_pw( pt )
1939             ENDIF
1940          ELSE
1941             CALL advec_s_up( pt )
1942          ENDIF
1943       ENDIF
1944
1945       CALL diffusion_s( pt,                                                   &
1946                         surf_def_h(0)%shf, surf_def_h(1)%shf,                 &
1947                         surf_def_h(2)%shf,                                    &
1948                         surf_lsm_h%shf,    surf_usm_h%shf,                    &
1949                         surf_def_v(0)%shf, surf_def_v(1)%shf,                 &
1950                         surf_def_v(2)%shf, surf_def_v(3)%shf,                 &
1951                         surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,                 &
1952                         surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,                 &
1953                         surf_usm_v(0)%shf, surf_usm_v(1)%shf,                 &
1954                         surf_usm_v(2)%shf, surf_usm_v(3)%shf )
1955
1956!
1957!--    Consideration of heat sources within the plant canopy
1958       IF ( plant_canopy  .AND.                                          &
1959            (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
1960          CALL pcm_tendency( 4 )
1961       ENDIF
1962
1963!
1964!--    Large scale advection
1965       IF ( large_scale_forcing )  THEN
1966          CALL ls_advec( simulated_time, 'pt' )
1967       ENDIF
1968
1969!
1970!--    Nudging
1971       IF ( nudging )  CALL nudge( simulated_time, 'pt' )
1972
1973!
1974!--    If required compute influence of large-scale subsidence/ascent
1975       IF ( large_scale_subsidence  .AND.                                      &
1976            .NOT. use_subsidence_tendencies )  THEN
1977          CALL subsidence( tend, pt, pt_init, 2 )
1978       ENDIF
1979
1980#if defined( __rrtmg )
1981!
1982!--    If required, add tendency due to radiative heating/cooling
1983       IF ( radiation  .AND.                                                   &
1984            simulated_time > skip_time_do_radiation )  THEN
1985            CALL radiation_tendency ( tend )
1986       ENDIF
1987#endif
1988
1989       CALL module_interface_actions( 'pt-tendency' )
1990
1991!
1992!--    Prognostic equation for potential temperature
1993       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1994       !$ACC PRESENT(pt, tend, tpt_m, wall_flags_0) &
1995       !$ACC PRESENT(pt_init, rdf_sc, ptdf_x, ptdf_y) &
1996       !$ACC PRESENT(tsc(3:5)) &
1997       !$ACC PRESENT(pt_p)
1998       DO  i = nxl, nxr
1999          DO  j = nys, nyn
2000             DO  k = nzb+1, nzt
2001                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +      &
2002                                                   tsc(3) * tpt_m(k,j,i) )     &
2003                                                 - tsc(5) *                    &
2004                                                   ( pt(k,j,i) - pt_init(k) ) *&
2005                                          ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )&
2006                                          )                                    &
2007                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2008                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2009                                          )
2010             ENDDO
2011          ENDDO
2012       ENDDO
2013!
2014!--    Calculate tendencies for the next Runge-Kutta step
2015       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2016          IF ( intermediate_timestep_count == 1 )  THEN
2017             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
2018             !$ACC PRESENT(tend, tpt_m)
2019             DO  i = nxl, nxr
2020                DO  j = nys, nyn
2021                   DO  k = nzb+1, nzt
2022                      tpt_m(k,j,i) = tend(k,j,i)
2023                   ENDDO
2024                ENDDO
2025             ENDDO
2026          ELSEIF ( intermediate_timestep_count < &
2027                   intermediate_timestep_count_max )  THEN
2028             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
2029             !$ACC PRESENT(tend, tpt_m)
2030             DO  i = nxl, nxr
2031                DO  j = nys, nyn
2032                   DO  k = nzb+1, nzt
2033                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
2034                                        5.3125_wp * tpt_m(k,j,i)
2035                   ENDDO
2036                ENDDO
2037             ENDDO
2038          ENDIF
2039       ENDIF
2040
2041       CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
2042
2043    ENDIF
2044
2045!
2046!-- If required, compute prognostic equation for total water content
2047    IF ( humidity )  THEN
2048
2049       CALL cpu_log( log_point(29), 'q-equation', 'start' )
2050
2051!
2052!--    Scalar/q-tendency terms with communication
2053       sbt = tsc(2)
2054       IF ( scalar_advec == 'bc-scheme' )  THEN
2055
2056          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2057!
2058!--          Bott-Chlond scheme always uses Euler time step. Thus:
2059             sbt = 1.0_wp
2060          ENDIF
2061          tend = 0.0_wp
2062          CALL advec_s_bc( q, 'q' )
2063
2064       ENDIF
2065
2066!
2067!--    Scalar/q-tendency terms with no communication
2068       IF ( scalar_advec /= 'bc-scheme' )  THEN
2069          tend = 0.0_wp
2070          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2071             IF ( ws_scheme_sca )  THEN
2072                CALL advec_s_ws( q, 'q' )
2073             ELSE
2074                CALL advec_s_pw( q )
2075             ENDIF
2076          ELSE
2077             CALL advec_s_up( q )
2078          ENDIF
2079       ENDIF
2080
2081       CALL diffusion_s( q,                                                    &
2082                         surf_def_h(0)%qsws, surf_def_h(1)%qsws,               &
2083                         surf_def_h(2)%qsws,                                   &
2084                         surf_lsm_h%qsws,    surf_usm_h%qsws,                  &
2085                         surf_def_v(0)%qsws, surf_def_v(1)%qsws,               &
2086                         surf_def_v(2)%qsws, surf_def_v(3)%qsws,               &
2087                         surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,               &
2088                         surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,               &
2089                         surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,               &
2090                         surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
2091
2092!
2093!--    Sink or source of humidity due to canopy elements
2094       IF ( plant_canopy ) CALL pcm_tendency( 5 )
2095
2096!
2097!--    Large scale advection
2098       IF ( large_scale_forcing )  THEN
2099          CALL ls_advec( simulated_time, 'q' )
2100       ENDIF
2101
2102!
2103!--    Nudging
2104       IF ( nudging )  CALL nudge( simulated_time, 'q' )
2105
2106!
2107!--    If required compute influence of large-scale subsidence/ascent
2108       IF ( large_scale_subsidence  .AND.                                      &
2109            .NOT. use_subsidence_tendencies )  THEN
2110         CALL subsidence( tend, q, q_init, 3 )
2111       ENDIF
2112
2113       CALL module_interface_actions( 'q-tendency' )
2114
2115!
2116!--    Prognostic equation for total water content
2117       DO  i = nxl, nxr
2118          DO  j = nys, nyn
2119             DO  k = nzb+1, nzt
2120                q_p(k,j,i) = q(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
2121                                                 tsc(3) * tq_m(k,j,i) )        &
2122                                               - tsc(5) * rdf_sc(k) *          &
2123                                                      ( q(k,j,i) - q_init(k) ) &
2124                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
2125                                                   BTEST( wall_flags_0(k,j,i), 0 ) &
2126                                                 )
2127                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
2128             ENDDO
2129          ENDDO
2130       ENDDO
2131
2132!
2133!--    Calculate tendencies for the next Runge-Kutta step
2134       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2135          IF ( intermediate_timestep_count == 1 )  THEN
2136             DO  i = nxl, nxr
2137                DO  j = nys, nyn
2138                   DO  k = nzb+1, nzt
2139                      tq_m(k,j,i) = tend(k,j,i)
2140                   ENDDO
2141                ENDDO
2142             ENDDO
2143          ELSEIF ( intermediate_timestep_count < &
2144                   intermediate_timestep_count_max )  THEN
2145             DO  i = nxl, nxr
2146                DO  j = nys, nyn
2147                   DO  k = nzb+1, nzt
2148                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
2149                                     + 5.3125_wp * tq_m(k,j,i)
2150                   ENDDO
2151                ENDDO
2152             ENDDO
2153          ENDIF
2154       ENDIF
2155
2156       CALL cpu_log( log_point(29), 'q-equation', 'stop' )
2157
2158!
2159!--    If required, calculate prognostic equations for cloud water content
2160!--    and cloud drop concentration
2161       IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
2162
2163          CALL cpu_log( log_point(67), 'qc-equation', 'start' )
2164
2165!
2166!--       Calculate prognostic equation for cloud water content
2167          sbt = tsc(2)
2168          IF ( scalar_advec == 'bc-scheme' )  THEN
2169
2170             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2171!
2172!--             Bott-Chlond scheme always uses Euler time step. Thus:
2173                sbt = 1.0_wp
2174             ENDIF
2175             tend = 0.0_wp
2176             CALL advec_s_bc( qc, 'qc' )
2177
2178          ENDIF
2179
2180!
2181!--       qc-tendency terms with no communication
2182          IF ( scalar_advec /= 'bc-scheme' )  THEN
2183             tend = 0.0_wp
2184             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2185                IF ( ws_scheme_sca )  THEN
2186                   CALL advec_s_ws( qc, 'qc' )
2187                ELSE
2188                   CALL advec_s_pw( qc )
2189                ENDIF
2190             ELSE
2191                CALL advec_s_up( qc )
2192             ENDIF
2193          ENDIF
2194
2195          CALL diffusion_s( qc,                                                &
2196                            surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,          &
2197                            surf_def_h(2)%qcsws,                               &
2198                            surf_lsm_h%qcsws,    surf_usm_h%qcsws,             &
2199                            surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,          &
2200                            surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,          &
2201                            surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,          &
2202                            surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,          &
2203                            surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,          &
2204                            surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
2205
2206!
2207!--       Prognostic equation for cloud water content
2208          DO  i = nxl, nxr
2209             DO  j = nys, nyn
2210                DO  k = nzb+1, nzt
2211                   qc_p(k,j,i) = qc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
2212                                                      tsc(3) * tqc_m(k,j,i) )  &
2213                                                    - tsc(5) * rdf_sc(k) *     &
2214                                                               qc(k,j,i)       &
2215                                             )                                 &
2216                                    * MERGE( 1.0_wp, 0.0_wp,                   &
2217                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2218                                          )
2219                   IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
2220                ENDDO
2221             ENDDO
2222          ENDDO
2223
2224!
2225!--       Calculate tendencies for the next Runge-Kutta step
2226          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2227             IF ( intermediate_timestep_count == 1 )  THEN
2228                DO  i = nxl, nxr
2229                   DO  j = nys, nyn
2230                      DO  k = nzb+1, nzt
2231                         tqc_m(k,j,i) = tend(k,j,i)
2232                      ENDDO
2233                   ENDDO
2234                ENDDO
2235             ELSEIF ( intermediate_timestep_count < &
2236                      intermediate_timestep_count_max )  THEN
2237                DO  i = nxl, nxr
2238                   DO  j = nys, nyn
2239                      DO  k = nzb+1, nzt
2240                         tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i)             &
2241                                         + 5.3125_wp * tqc_m(k,j,i)
2242                      ENDDO
2243                   ENDDO
2244                ENDDO
2245             ENDIF
2246          ENDIF
2247
2248          CALL cpu_log( log_point(67), 'qc-equation', 'stop' )
2249
2250          CALL cpu_log( log_point(68), 'nc-equation', 'start' )
2251!
2252!--       Calculate prognostic equation for cloud drop concentration
2253          sbt = tsc(2)
2254          IF ( scalar_advec == 'bc-scheme' )  THEN
2255
2256             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2257!
2258!--             Bott-Chlond scheme always uses Euler time step. Thus:
2259                sbt = 1.0_wp
2260             ENDIF
2261             tend = 0.0_wp
2262             CALL advec_s_bc( nc, 'nc' )
2263
2264          ENDIF
2265
2266!
2267!--       nc-tendency terms with no communication
2268          IF ( scalar_advec /= 'bc-scheme' )  THEN
2269             tend = 0.0_wp
2270             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2271                IF ( ws_scheme_sca )  THEN
2272                   CALL advec_s_ws( nc, 'nc' )
2273                ELSE
2274                   CALL advec_s_pw( nc )
2275                ENDIF
2276             ELSE
2277                CALL advec_s_up( nc )
2278             ENDIF
2279          ENDIF
2280
2281          CALL diffusion_s( nc,                                                &
2282                            surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,          &
2283                            surf_def_h(2)%ncsws,                               &
2284                            surf_lsm_h%ncsws,    surf_usm_h%ncsws,             & 
2285                            surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,          &
2286                            surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,          &
2287                            surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,          &
2288                            surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,          &
2289                            surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,          &
2290                            surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
2291
2292!
2293!--       Prognostic equation for cloud drop concentration
2294          DO  i = nxl, nxr
2295             DO  j = nys, nyn
2296                DO  k = nzb+1, nzt
2297                   nc_p(k,j,i) = nc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
2298                                                      tsc(3) * tnc_m(k,j,i) )  &
2299                                                    - tsc(5) * rdf_sc(k) *     &
2300                                                               nc(k,j,i)       &
2301                                             )                                 &
2302                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2303                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2304                                          )
2305                   IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
2306                ENDDO
2307             ENDDO
2308          ENDDO
2309
2310!
2311!--       Calculate tendencies for the next Runge-Kutta step
2312          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2313             IF ( intermediate_timestep_count == 1 )  THEN
2314                DO  i = nxl, nxr
2315                   DO  j = nys, nyn
2316                      DO  k = nzb+1, nzt
2317                         tnc_m(k,j,i) = tend(k,j,i)
2318                      ENDDO
2319                   ENDDO
2320                ENDDO
2321             ELSEIF ( intermediate_timestep_count < &
2322                      intermediate_timestep_count_max )  THEN
2323                DO  i = nxl, nxr
2324                   DO  j = nys, nyn
2325                      DO  k = nzb+1, nzt
2326                         tnc_m(k,j,i) =  -9.5625_wp * tend(k,j,i)             &
2327                                         + 5.3125_wp * tnc_m(k,j,i)
2328                      ENDDO
2329                   ENDDO
2330                ENDDO
2331             ENDIF
2332          ENDIF
2333
2334          CALL cpu_log( log_point(68), 'nc-equation', 'stop' )
2335
2336       ENDIF
2337!
2338!--    If required, calculate prognostic equations for rain water content
2339!--    and rain drop concentration
2340       IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
2341
2342          CALL cpu_log( log_point(52), 'qr-equation', 'start' )
2343
2344!
2345!--       Calculate prognostic equation for rain water content
2346          sbt = tsc(2)
2347          IF ( scalar_advec == 'bc-scheme' )  THEN
2348
2349             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2350!
2351!--             Bott-Chlond scheme always uses Euler time step. Thus:
2352                sbt = 1.0_wp
2353             ENDIF
2354             tend = 0.0_wp
2355             CALL advec_s_bc( qr, 'qr' )
2356
2357          ENDIF
2358
2359!
2360!--       qr-tendency terms with no communication
2361          IF ( scalar_advec /= 'bc-scheme' )  THEN
2362             tend = 0.0_wp
2363             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2364                IF ( ws_scheme_sca )  THEN
2365                   CALL advec_s_ws( qr, 'qr' )
2366                ELSE
2367                   CALL advec_s_pw( qr )
2368                ENDIF
2369             ELSE
2370                CALL advec_s_up( qr )
2371             ENDIF
2372          ENDIF
2373
2374          CALL diffusion_s( qr,                                                &
2375                            surf_def_h(0)%qrsws, surf_def_h(1)%qrsws,          &
2376                            surf_def_h(2)%qrsws,                               &
2377                            surf_lsm_h%qrsws,    surf_usm_h%qrsws,             &
2378                            surf_def_v(0)%qrsws, surf_def_v(1)%qrsws,          &
2379                            surf_def_v(2)%qrsws, surf_def_v(3)%qrsws,          &
2380                            surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws,          &
2381                            surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws,          &
2382                            surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws,          &
2383                            surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws )
2384
2385!
2386!--       Prognostic equation for rain water content
2387          DO  i = nxl, nxr
2388             DO  j = nys, nyn
2389                DO  k = nzb+1, nzt
2390                   qr_p(k,j,i) = qr(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
2391                                                      tsc(3) * tqr_m(k,j,i) )  &
2392                                                    - tsc(5) * rdf_sc(k) *     &
2393                                                               qr(k,j,i)       &
2394                                             )                                 &
2395                                    * MERGE( 1.0_wp, 0.0_wp,                   &
2396                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2397                                          )
2398                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
2399                ENDDO
2400             ENDDO
2401          ENDDO
2402
2403!
2404!--       Calculate tendencies for the next Runge-Kutta step
2405          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2406             IF ( intermediate_timestep_count == 1 )  THEN
2407                DO  i = nxl, nxr
2408                   DO  j = nys, nyn
2409                      DO  k = nzb+1, nzt
2410                         tqr_m(k,j,i) = tend(k,j,i)
2411                      ENDDO
2412                   ENDDO
2413                ENDDO
2414             ELSEIF ( intermediate_timestep_count < &
2415                      intermediate_timestep_count_max )  THEN
2416                DO  i = nxl, nxr
2417                   DO  j = nys, nyn
2418                      DO  k = nzb+1, nzt
2419                         tqr_m(k,j,i) =   -9.5625_wp * tend(k,j,i)             &
2420                                         + 5.3125_wp * tqr_m(k,j,i)
2421                      ENDDO
2422                   ENDDO
2423                ENDDO
2424             ENDIF
2425          ENDIF
2426
2427          CALL cpu_log( log_point(52), 'qr-equation', 'stop' )
2428          CALL cpu_log( log_point(53), 'nr-equation', 'start' )
2429
2430!
2431!--       Calculate prognostic equation for rain drop concentration
2432          sbt = tsc(2)
2433          IF ( scalar_advec == 'bc-scheme' )  THEN
2434
2435             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2436!
2437!--             Bott-Chlond scheme always uses Euler time step. Thus:
2438                sbt = 1.0_wp
2439             ENDIF
2440             tend = 0.0_wp
2441             CALL advec_s_bc( nr, 'nr' )
2442
2443          ENDIF
2444
2445!
2446!--       nr-tendency terms with no communication
2447          IF ( scalar_advec /= 'bc-scheme' )  THEN
2448             tend = 0.0_wp
2449             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2450                IF ( ws_scheme_sca )  THEN
2451                   CALL advec_s_ws( nr, 'nr' )
2452                ELSE
2453                   CALL advec_s_pw( nr )
2454                ENDIF
2455             ELSE
2456                CALL advec_s_up( nr )
2457             ENDIF
2458          ENDIF
2459
2460          CALL diffusion_s( nr,                                                &
2461                            surf_def_h(0)%nrsws, surf_def_h(1)%nrsws,          &
2462                            surf_def_h(2)%nrsws,                               &
2463                            surf_lsm_h%nrsws,    surf_usm_h%nrsws,             & 
2464                            surf_def_v(0)%nrsws, surf_def_v(1)%nrsws,          &
2465                            surf_def_v(2)%nrsws, surf_def_v(3)%nrsws,          &
2466                            surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws,          &
2467                            surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws,          &
2468                            surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws,          &
2469                            surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws )
2470
2471!
2472!--       Prognostic equation for rain drop concentration
2473          DO  i = nxl, nxr
2474             DO  j = nys, nyn
2475                DO  k = nzb+1, nzt
2476                   nr_p(k,j,i) = nr(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
2477                                                      tsc(3) * tnr_m(k,j,i) )  &
2478                                                    - tsc(5) * rdf_sc(k) *     &
2479                                                               nr(k,j,i)       &
2480                                             )                                 &
2481                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2482                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2483                                          )
2484                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
2485                ENDDO
2486             ENDDO
2487          ENDDO
2488
2489!
2490!--       Calculate tendencies for the next Runge-Kutta step
2491          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2492             IF ( intermediate_timestep_count == 1 )  THEN
2493                DO  i = nxl, nxr
2494                   DO  j = nys, nyn
2495                      DO  k = nzb+1, nzt
2496                         tnr_m(k,j,i) = tend(k,j,i)
2497                      ENDDO
2498                   ENDDO
2499                ENDDO
2500             ELSEIF ( intermediate_timestep_count < &
2501                      intermediate_timestep_count_max )  THEN
2502                DO  i = nxl, nxr
2503                   DO  j = nys, nyn
2504                      DO  k = nzb+1, nzt
2505                         tnr_m(k,j,i) =  -9.5625_wp * tend(k,j,i)             &
2506                                         + 5.3125_wp * tnr_m(k,j,i)
2507                      ENDDO
2508                   ENDDO
2509                ENDDO
2510             ENDIF
2511          ENDIF
2512
2513          CALL cpu_log( log_point(53), 'nr-equation', 'stop' )
2514
2515       ENDIF
2516
2517    ENDIF
2518!
2519!-- If required, compute prognostic equation for scalar
2520    IF ( passive_scalar )  THEN
2521
2522       CALL cpu_log( log_point(66), 's-equation', 'start' )
2523
2524!
2525!--    Scalar/q-tendency terms with communication
2526       sbt = tsc(2)
2527       IF ( scalar_advec == 'bc-scheme' )  THEN
2528
2529          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2530!
2531!--          Bott-Chlond scheme always uses Euler time step. Thus:
2532             sbt = 1.0_wp
2533          ENDIF
2534          tend = 0.0_wp
2535          CALL advec_s_bc( s, 's' )
2536
2537       ENDIF
2538
2539!
2540!--    Scalar/q-tendency terms with no communication
2541       IF ( scalar_advec /= 'bc-scheme' )  THEN
2542          tend = 0.0_wp
2543          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2544             IF ( ws_scheme_sca )  THEN
2545                CALL advec_s_ws( s, 's' )
2546             ELSE
2547                CALL advec_s_pw( s )
2548             ENDIF
2549          ELSE
2550             CALL advec_s_up( s )
2551          ENDIF
2552       ENDIF
2553
2554       CALL diffusion_s( s,                                                    &
2555                         surf_def_h(0)%ssws, surf_def_h(1)%ssws,               &
2556                         surf_def_h(2)%ssws,                                   &
2557                         surf_lsm_h%ssws,    surf_usm_h%ssws,                  &
2558                         surf_def_v(0)%ssws, surf_def_v(1)%ssws,               &
2559                         surf_def_v(2)%ssws, surf_def_v(3)%ssws,               &
2560                         surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,               &
2561                         surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,               &
2562                         surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,               &
2563                         surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
2564
2565!
2566!--    Sink or source of humidity due to canopy elements
2567       IF ( plant_canopy ) CALL pcm_tendency( 7 )
2568
2569!
2570!--    Large scale advection. Not implemented for scalars so far.
2571!        IF ( large_scale_forcing )  THEN
2572!           CALL ls_advec( simulated_time, 'q' )
2573!        ENDIF
2574
2575!
2576!--    Nudging. Not implemented for scalars so far.
2577!        IF ( nudging )  CALL nudge( simulated_time, 'q' )
2578
2579!
2580!--    If required compute influence of large-scale subsidence/ascent.
2581!--    Not implemented for scalars so far.
2582       IF ( large_scale_subsidence  .AND.                                      &
2583            .NOT. use_subsidence_tendencies  .AND.                             &
2584            .NOT. large_scale_forcing )  THEN
2585         CALL subsidence( tend, s, s_init, 3 )
2586       ENDIF
2587
2588       CALL module_interface_actions( 's-tendency' )
2589
2590!
2591!--    Prognostic equation for total water content
2592       DO  i = nxl, nxr
2593          DO  j = nys, nyn
2594             DO  k = nzb+1, nzt
2595                s_p(k,j,i) = s(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
2596                                                 tsc(3) * ts_m(k,j,i) )        &
2597                                               - tsc(5) * rdf_sc(k) *          &
2598                                                    ( s(k,j,i) - s_init(k) )   &
2599                                        )                                      &
2600                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2601                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2602                                          )
2603                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
2604             ENDDO
2605          ENDDO
2606       ENDDO
2607
2608!
2609!--    Calculate tendencies for the next Runge-Kutta step
2610       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2611          IF ( intermediate_timestep_count == 1 )  THEN
2612             DO  i = nxl, nxr
2613                DO  j = nys, nyn
2614                   DO  k = nzb+1, nzt
2615                      ts_m(k,j,i) = tend(k,j,i)
2616                   ENDDO
2617                ENDDO
2618             ENDDO
2619          ELSEIF ( intermediate_timestep_count < &
2620                   intermediate_timestep_count_max )  THEN
2621             DO  i = nxl, nxr
2622                DO  j = nys, nyn
2623                   DO  k = nzb+1, nzt
2624                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
2625                                     + 5.3125_wp * ts_m(k,j,i)
2626                   ENDDO
2627                ENDDO
2628             ENDDO
2629          ENDIF
2630       ENDIF
2631
2632       CALL cpu_log( log_point(66), 's-equation', 'stop' )
2633
2634    ENDIF
2635
2636!
2637!-- Calculate prognostic equations for turbulence closure
2638    CALL tcm_prognostic_equations()
2639
2640!
2641!-- Calculate prognostic equation for chemical quantites
2642    IF ( air_chemistry )  THEN
2643       CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'start' )
2644!
2645!--    Loop over chemical species
2646       DO  lsp = 1, nvar                         
2647          CALL chem_prognostic_equations( chem_species(lsp)%conc_p,            &
2648                                          chem_species(lsp)%conc,              &
2649                                          chem_species(lsp)%tconc_m,           &
2650                                          chem_species(lsp)%conc_pr_init,      &
2651                                          lsp )
2652       ENDDO
2653
2654       CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'stop' )             
2655    ENDIF   ! Chemicals equations
2656       
2657    IF ( salsa )  THEN
2658       CALL cpu_log( log_point_s(92), 'salsa advec+diff+prog ', 'start' )
2659!
2660!--          Loop over aerosol size bins: number and mass bins
2661       IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
2662       
2663          DO  b = 1, nbins               
2664             sums_salsa_ws_l = aerosol_number(b)%sums_ws_l
2665             CALL salsa_tendency( 'aerosol_number', aerosol_number(b)%conc_p,  &
2666                                   aerosol_number(b)%conc,                     &
2667                                   aerosol_number(b)%tconc_m,                  &
2668                                   b, b, aerosol_number(b)%init )
2669             aerosol_number(b)%sums_ws_l = sums_salsa_ws_l
2670             DO  c = 1, ncc_tot
2671                sums_salsa_ws_l = aerosol_mass((c-1)*nbins+b)%sums_ws_l 
2672                CALL salsa_tendency( 'aerosol_mass',                           &
2673                                      aerosol_mass((c-1)*nbins+b)%conc_p,      &
2674                                      aerosol_mass((c-1)*nbins+b)%conc,        &
2675                                      aerosol_mass((c-1)*nbins+b)%tconc_m,     &
2676                                      b, c, aerosol_mass((c-1)*nbins+b)%init )
2677                aerosol_mass((c-1)*nbins+b)%sums_ws_l = sums_salsa_ws_l
2678             ENDDO
2679          ENDDO
2680          IF ( .NOT. salsa_gases_from_chem )  THEN
2681             DO  g = 1, ngast
2682                sums_salsa_ws_l = salsa_gas(g)%sums_ws_l
2683                CALL salsa_tendency( 'salsa_gas', salsa_gas(g)%conc_p,         &
2684                                      salsa_gas(g)%conc, salsa_gas(g)%tconc_m, &
2685                                      g, g, salsa_gas(g)%init )
2686                salsa_gas(g)%sums_ws_l = sums_salsa_ws_l
2687             ENDDO
2688          ENDIF
2689       
2690       ENDIF
2691       
2692       CALL cpu_log( log_point_s(92), 'salsa advec+diff+prog ', 'stop' )
2693    ENDIF 
2694
2695!
2696!-- Calculate prognostic equations for the ocean
2697    IF ( ocean_mode )  CALL ocean_prognostic_equations()
2698
2699 END SUBROUTINE prognostic_equations_vector
2700
2701
2702 END MODULE prognostic_equations_mod
Note: See TracBrowser for help on using the repository browser.