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

Last change on this file since 3804 was 3797, checked in by forkel, 5 years ago

Modifications for OpenMP version by Klaus Ketelsen

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