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

Last change on this file since 2232 was 2232, checked in by suehring, 4 years ago

Adjustments according new topography and surface-modelling concept implemented

  • Property svn:keywords set to Id
File size: 81.9 KB
Line 
1!> @file prognostic_equations.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
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-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22! Adjutst to new surface-type structure. Remove call for usm_wall_heat_flux,
23! which is realized directly in diffusion_s now.
24!
25! Former revisions:
26! -----------------
27! $Id: prognostic_equations.f90 2232 2017-05-30 17:47:52Z suehring $
28!
29! 2192 2017-03-22 04:14:10Z raasch
30! Bugfix for misplaced and missing openMP directives from r2155
31!
32! 2155 2017-02-21 09:57:40Z hoffmann
33! Bugfix in the calculation of microphysical quantities on ghost points.
34!
35! 2118 2017-01-17 16:38:49Z raasch
36! OpenACC version of subroutine removed
37!
38! 2031 2016-10-21 15:11:58Z knoop
39! renamed variable rho to rho_ocean
40!
41! 2011 2016-09-19 17:29:57Z kanani
42! Flag urban_surface is now defined in module control_parameters.
43!
44! 2007 2016-08-24 15:47:17Z kanani
45! Added pt tendency calculation based on energy balance at urban surfaces
46! (new urban surface model)
47!
48! 2000 2016-08-20 18:09:15Z knoop
49! Forced header and separation lines into 80 columns
50!
51! 1976 2016-07-27 13:28:04Z maronga
52! Simplied calls to radiation model
53!
54! 1960 2016-07-12 16:34:24Z suehring
55! Separate humidity and passive scalar
56!
57! 1914 2016-05-26 14:44:07Z witha
58! Added calls for wind turbine model
59!
60! 1873 2016-04-18 14:50:06Z maronga
61! Module renamed (removed _mod)
62!
63! 1850 2016-04-08 13:29:27Z maronga
64! Module renamed
65!
66! 1826 2016-04-07 12:01:39Z maronga
67! Renamed canopy model calls.
68!
69! 1822 2016-04-07 07:49:42Z hoffmann
70! Kessler microphysics scheme moved to microphysics.
71!
72! 1757 2016-02-22 15:49:32Z maronga
73!
74! 1691 2015-10-26 16:17:44Z maronga
75! Added optional model spin-up without radiation / land surface model calls.
76! Formatting corrections.
77!
78! 1682 2015-10-07 23:56:08Z knoop
79! Code annotations made doxygen readable
80!
81! 1585 2015-04-30 07:05:52Z maronga
82! Added call for temperature tendency calculation due to radiative flux divergence
83!
84! 1517 2015-01-07 19:12:25Z hoffmann
85! advec_s_bc_mod addded, since advec_s_bc is now a module
86!
87! 1496 2014-12-02 17:25:50Z maronga
88! Renamed "radiation" -> "cloud_top_radiation"
89!
90! 1484 2014-10-21 10:53:05Z kanani
91! Changes due to new module structure of the plant canopy model:
92!   parameters cthf and plant_canopy moved to module plant_canopy_model_mod.
93! Removed double-listing of use_upstream_for_tke in ONLY-list of module
94! control_parameters
95!
96! 1409 2014-05-23 12:11:32Z suehring
97! Bugfix: i_omp_start changed for advec_u_ws at left inflow and outflow boundary.
98! This ensures that left-hand side fluxes are also calculated for nxl in that
99! case, even though the solution at nxl is overwritten in boundary_conds()
100!
101! 1398 2014-05-07 11:15:00Z heinze
102! Rayleigh-damping for horizontal velocity components changed: instead of damping
103! against ug and vg, damping against u_init and v_init is used to allow for a
104! homogenized treatment in case of nudging
105!
106! 1380 2014-04-28 12:40:45Z heinze
107! Change order of calls for scalar prognostic quantities:
108! ls_advec -> nudging -> subsidence since initial profiles
109!
110! 1374 2014-04-25 12:55:07Z raasch
111! missing variables added to ONLY lists
112!
113! 1365 2014-04-22 15:03:56Z boeske
114! Calls of ls_advec for large scale advection added,
115! subroutine subsidence is only called if use_subsidence_tendencies = .F.,
116! new argument ls_index added to the calls of subsidence
117! +ls_index
118!
119! 1361 2014-04-16 15:17:48Z hoffmann
120! Two-moment microphysics moved to the start of prognostic equations. This makes
121! the 3d arrays for tend_q, tend_qr, tend_pt and tend_pt redundant.
122! Additionally, it is allowed to call the microphysics just once during the time
123! step (not at each sub-time step).
124!
125! Two-moment cloud physics added for vector and accelerator optimization.
126!
127! 1353 2014-04-08 15:21:23Z heinze
128! REAL constants provided with KIND-attribute
129!
130! 1337 2014-03-25 15:11:48Z heinze
131! Bugfix: REAL constants provided with KIND-attribute
132!
133! 1332 2014-03-25 11:59:43Z suehring
134! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke
135!
136! 1330 2014-03-24 17:29:32Z suehring
137! In case of SGS-particle velocity advection of TKE is also allowed with
138! dissipative 5th-order scheme.
139!
140! 1320 2014-03-20 08:40:49Z raasch
141! ONLY-attribute added to USE-statements,
142! kind-parameters added to all INTEGER and REAL declaration statements,
143! kinds are defined in new module kinds,
144! old module precision_kind is removed,
145! revision history before 2012 removed,
146! comment fields (!:) to be used for variable explanations added to
147! all variable declaration statements
148!
149! 1318 2014-03-17 13:35:16Z raasch
150! module interfaces removed
151!
152! 1257 2013-11-08 15:18:40Z raasch
153! openacc loop vector clauses removed, independent clauses added
154!
155! 1246 2013-11-01 08:59:45Z heinze
156! enable nudging also for accelerator version
157!
158! 1241 2013-10-30 11:36:58Z heinze
159! usage of nudging enabled (so far not implemented for accelerator version)
160!
161! 1179 2013-06-14 05:57:58Z raasch
162! two arguments removed from routine buoyancy, ref_state updated on device
163!
164! 1128 2013-04-12 06:19:32Z raasch
165! those parts requiring global communication moved to time_integration,
166! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
167! j_north
168!
169! 1115 2013-03-26 18:16:16Z hoffmann
170! optimized cloud physics: calculation of microphysical tendencies transfered
171! to microphysics.f90; qr and nr are only calculated if precipitation is required
172!
173! 1111 2013-03-08 23:54:10Z raasch
174! update directives for prognostic quantities removed
175!
176! 1106 2013-03-04 05:31:38Z raasch
177! small changes in code formatting
178!
179! 1092 2013-02-02 11:24:22Z raasch
180! unused variables removed
181!
182! 1053 2012-11-13 17:11:03Z hoffmann
183! implementation of two new prognostic equations for rain drop concentration (nr)
184! and rain water content (qr)
185!
186! currently, only available for cache loop optimization
187!
188! 1036 2012-10-22 13:43:42Z raasch
189! code put under GPL (PALM 3.9)
190!
191! 1019 2012-09-28 06:46:45Z raasch
192! non-optimized version of prognostic_equations removed
193!
194! 1015 2012-09-27 09:23:24Z raasch
195! new branch prognostic_equations_acc
196! OpenACC statements added + code changes required for GPU optimization
197!
198! 1001 2012-09-13 14:08:46Z raasch
199! all actions concerning leapfrog- and upstream-spline-scheme removed
200!
201! 978 2012-08-09 08:28:32Z fricke
202! km_damp_x and km_damp_y removed in calls of diffusion_u and diffusion_v
203! add ptdf_x, ptdf_y for damping the potential temperature at the inflow
204! boundary in case of non-cyclic lateral boundaries
205! Bugfix: first thread index changes for WS-scheme at the inflow
206!
207! 940 2012-07-09 14:31:00Z raasch
208! temperature equation can be switched off
209!
210! Revision 1.1  2000/04/13 14:56:27  schroeter
211! Initial revision
212!
213!
214! Description:
215! ------------
216!> Solving the prognostic equations.
217!------------------------------------------------------------------------------!
218 MODULE prognostic_equations_mod
219
220
221
222    USE arrays_3d,                                                             &
223        ONLY:  diss_l_e, diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr,            &
224               diss_l_s, diss_l_sa, diss_s_e, diss_s_nr, diss_s_pt, diss_s_q,  &
225               diss_s_qr, diss_s_s, diss_s_sa, e, e_p, flux_s_e, flux_s_nr,    &
226               flux_s_pt, flux_s_q, flux_s_qr, flux_s_s, flux_s_sa, flux_l_e,  &
227               flux_l_nr, flux_l_pt, flux_l_q, flux_l_qr, flux_l_s, flux_l_sa, &
228               nr, nr_p, pt, ptdf_x, ptdf_y, pt_init, pt_p,                    &
229               prho, q, q_init, q_p, qr, qr_p, rdf,                            &
230               rdf_sc, ref_state, rho_ocean, s, s_init, s_p, sa, sa_init, sa_p,&
231               tend, te_m, tnr_m, tpt_m, tq_m, tqr_m, ts_m, tsa_m, tu_m, tv_m, &
232               tw_m, u, ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
233
234    USE control_parameters,                                                    &
235        ONLY:  call_microphysics_at_all_substeps, cloud_physics,               &
236               cloud_top_radiation, constant_diffusion, dp_external,           &
237               dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d, humidity,       &
238               inflow_l, intermediate_timestep_count,                          &
239               intermediate_timestep_count_max, large_scale_forcing,           &
240               large_scale_subsidence, microphysics_seifert,                   &
241               microphysics_sat_adjust, neutral, nudging, ocean, outflow_l,    &
242               outflow_s, passive_scalar, prho_reference, prho_reference,      &
243               prho_reference, pt_reference, pt_reference, pt_reference,       &
244               scalar_advec, scalar_advec, simulated_time, sloping_surface,    &
245               timestep_scheme, tsc, use_subsidence_tendencies,                &
246               use_upstream_for_tke, ws_scheme_mom, ws_scheme_sca
247
248    USE cpulog,                                                                &
249        ONLY:  cpu_log, log_point
250
251    USE eqn_state_seawater_mod,                                                &
252        ONLY:  eqn_state_seawater
253
254    USE indices,                                                               &
255        ONLY:  nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv, nzb,    &
256               nzt, wall_flags_0
257
258    USE advec_ws,                                                              &
259        ONLY:  advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws
260
261    USE advec_s_bc_mod,                                                        &
262        ONLY:  advec_s_bc
263
264    USE advec_s_pw_mod,                                                        &
265        ONLY:  advec_s_pw
266
267    USE advec_s_up_mod,                                                        &
268        ONLY:  advec_s_up
269
270    USE advec_u_pw_mod,                                                        &
271        ONLY:  advec_u_pw
272
273    USE advec_u_up_mod,                                                        &
274        ONLY:  advec_u_up
275
276    USE advec_v_pw_mod,                                                        &
277        ONLY:  advec_v_pw
278
279    USE advec_v_up_mod,                                                        &
280        ONLY:  advec_v_up
281
282    USE advec_w_pw_mod,                                                        &
283        ONLY:  advec_w_pw
284
285    USE advec_w_up_mod,                                                        &
286        ONLY:  advec_w_up
287
288    USE buoyancy_mod,                                                          &
289        ONLY:  buoyancy
290
291    USE calc_radiation_mod,                                                    &
292        ONLY:  calc_radiation
293
294    USE coriolis_mod,                                                          &
295        ONLY:  coriolis
296
297    USE diffusion_e_mod,                                                       &
298        ONLY:  diffusion_e
299
300    USE diffusion_s_mod,                                                       &
301        ONLY:  diffusion_s
302
303    USE diffusion_u_mod,                                                       &
304        ONLY:  diffusion_u
305
306    USE diffusion_v_mod,                                                       &
307        ONLY:  diffusion_v
308
309    USE diffusion_w_mod,                                                       &
310        ONLY:  diffusion_w
311
312    USE kinds
313
314    USE ls_forcing_mod,                                                        &
315        ONLY:  ls_advec
316
317    USE microphysics_mod,                                                      &
318        ONLY:  microphysics_control
319
320    USE nudge_mod,                                                             &
321        ONLY:  nudge
322
323    USE plant_canopy_model_mod,                                                &
324        ONLY:  cthf, plant_canopy, pcm_tendency
325
326    USE production_e_mod,                                                      &
327        ONLY:  production_e
328
329    USE radiation_model_mod,                                                   &
330        ONLY:  radiation, radiation_tendency,                                  &
331               skip_time_do_radiation
332
333    USE statistics,                                                            &
334        ONLY:  hom
335
336    USE subsidence_mod,                                                        &
337        ONLY:  subsidence
338
339    USE user_actions_mod,                                                      &
340        ONLY:  user_actions
341
342    USE surface_mod,                                                           &
343        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
344                surf_usm_v
345
346    USE wind_turbine_model_mod,                                                &
347        ONLY:  wind_turbine, wtm_tendencies
348
349
350    PRIVATE
351    PUBLIC prognostic_equations_cache, prognostic_equations_vector
352
353    INTERFACE prognostic_equations_cache
354       MODULE PROCEDURE prognostic_equations_cache
355    END INTERFACE prognostic_equations_cache
356
357    INTERFACE prognostic_equations_vector
358       MODULE PROCEDURE prognostic_equations_vector
359    END INTERFACE prognostic_equations_vector
360
361
362 CONTAINS
363
364
365!------------------------------------------------------------------------------!
366! Description:
367! ------------
368!> Version with one optimized loop over all equations. It is only allowed to
369!> be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
370!>
371!> Here the calls of most subroutines are embedded in two DO loops over i and j,
372!> so communication between CPUs is not allowed (does not make sense) within
373!> these loops.
374!>
375!> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
376!------------------------------------------------------------------------------!
377
378 SUBROUTINE prognostic_equations_cache
379
380
381    IMPLICIT NONE
382
383    INTEGER(iwp) ::  i                   !<
384    INTEGER(iwp) ::  i_omp_start         !<
385    INTEGER(iwp) ::  j                   !<
386    INTEGER(iwp) ::  k                   !<
387    INTEGER(iwp) ::  omp_get_thread_num  !<
388    INTEGER(iwp) ::  tn = 0              !<
389
390    LOGICAL      ::  loop_start          !<
391
392
393!
394!-- Time measurement can only be performed for the whole set of equations
395    CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
396
397!
398!-- If required, calculate cloud microphysics
399    IF ( cloud_physics  .AND.  .NOT. microphysics_sat_adjust  .AND.            &
400         ( intermediate_timestep_count == 1  .OR.                              &
401           call_microphysics_at_all_substeps ) )                               &
402    THEN
403       !$OMP PARALLEL private (i,j)
404       !$OMP DO
405!$OMP PARALLEL PRIVATE (i, j)
406!$OMP DO
407       DO  i = nxlg, nxrg
408          DO  j = nysg, nyng
409             CALL microphysics_control( i, j )
410           ENDDO
411       ENDDO
412       !$OMP END PARALLEL
413!$OMP END PARALLEL
414    ENDIF
415
416!
417!-- Loop over all prognostic equations
418    !$OMP PARALLEL private (i,i_omp_start,j,k,loop_start,tn)
419
420    !$  tn = omp_get_thread_num()
421    loop_start = .TRUE.
422
423    !$OMP DO
424    DO  i = nxl, nxr
425
426!
427!--    Store the first loop index. It differs for each thread and is required
428!--    later in advec_ws
429       IF ( loop_start )  THEN
430          loop_start  = .FALSE.
431          i_omp_start = i
432       ENDIF
433
434       DO  j = nys, nyn
435!
436!--       Tendency terms for u-velocity component
437          IF ( .NOT. outflow_l  .OR.  i > nxl )  THEN
438
439             tend(:,j,i) = 0.0_wp
440             IF ( timestep_scheme(1:5) == 'runge' )  THEN
441                IF ( ws_scheme_mom )  THEN
442                   CALL advec_u_ws( i, j, i_omp_start, tn )
443                ELSE
444                   CALL advec_u_pw( i, j )
445                ENDIF
446             ELSE
447                CALL advec_u_up( i, j )
448             ENDIF
449             CALL diffusion_u( i, j )
450             CALL coriolis( i, j, 1 )
451             IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
452                CALL buoyancy( i, j, pt, 1 )
453             ENDIF
454
455!
456!--          Drag by plant canopy
457             IF ( plant_canopy )  CALL pcm_tendency( i, j, 1 )
458
459!
460!--          External pressure gradient
461             IF ( dp_external )  THEN
462                DO  k = dp_level_ind_b+1, nzt
463                   tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
464                ENDDO
465             ENDIF
466
467!
468!--          Nudging
469             IF ( nudging )  CALL nudge( i, j, simulated_time, 'u' )
470
471!
472!--          Forces by wind turbines
473             IF ( wind_turbine )  CALL wtm_tendencies( i, j, 1 )
474
475             CALL user_actions( i, j, 'u-tendency' )
476!
477!--          Prognostic equation for u-velocity component
478             DO  k = nzb+1, nzt
479
480                u_p(k,j,i) = u(k,j,i) + ( dt_3d *                               &
481                                            ( tsc(2) * tend(k,j,i) +            &
482                                              tsc(3) * tu_m(k,j,i) )            &
483                                            - tsc(5) * rdf(k)                   &
484                                                     * ( u(k,j,i) - u_init(k) ) &
485                                        ) * MERGE( 1.0_wp, 0.0_wp,              &
486                                                   BTEST( wall_flags_0(k,j,i), 1 )&
487                                                 )
488             ENDDO
489
490!
491!--          Calculate tendencies for the next Runge-Kutta step
492             IF ( timestep_scheme(1:5) == 'runge' )  THEN
493                IF ( intermediate_timestep_count == 1 )  THEN
494                   DO  k = nzb+1, nzt
495                      tu_m(k,j,i) = tend(k,j,i)
496                   ENDDO
497                ELSEIF ( intermediate_timestep_count < &
498                         intermediate_timestep_count_max )  THEN
499                   DO  k = nzb+1, nzt
500                      tu_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                &
501                                     + 5.3125_wp * tu_m(k,j,i)
502                   ENDDO
503                ENDIF
504             ENDIF
505
506          ENDIF
507
508!
509!--       Tendency terms for v-velocity component
510          IF ( .NOT. outflow_s  .OR.  j > nys )  THEN
511
512             tend(:,j,i) = 0.0_wp
513             IF ( timestep_scheme(1:5) == 'runge' )  THEN
514                IF ( ws_scheme_mom )  THEN
515                    CALL advec_v_ws( i, j, i_omp_start, tn )
516                ELSE
517                    CALL advec_v_pw( i, j )
518                ENDIF
519             ELSE
520                CALL advec_v_up( i, j )
521             ENDIF
522             CALL diffusion_v( i, j )
523             CALL coriolis( i, j, 2 )
524
525!
526!--          Drag by plant canopy
527             IF ( plant_canopy )  CALL pcm_tendency( i, j, 2 )
528
529!
530!--          External pressure gradient
531             IF ( dp_external )  THEN
532                DO  k = dp_level_ind_b+1, nzt
533                   tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
534                ENDDO
535             ENDIF
536
537!
538!--          Nudging
539             IF ( nudging )  CALL nudge( i, j, simulated_time, 'v' )
540
541!
542!--          Forces by wind turbines
543             IF ( wind_turbine )  CALL wtm_tendencies( i, j, 2 )
544
545             CALL user_actions( i, j, 'v-tendency' )
546!
547!--          Prognostic equation for v-velocity component
548             DO  k = nzb+1, nzt
549                v_p(k,j,i) = v(k,j,i) + ( dt_3d *                              &
550                                            ( tsc(2) * tend(k,j,i) +           &
551                                              tsc(3) * tv_m(k,j,i) )           &
552                                            - tsc(5) * rdf(k)                  &
553                                                     * ( v(k,j,i) - v_init(k) )&
554                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
555                                                   BTEST( wall_flags_0(k,j,i), 2 )&
556                                                 )
557             ENDDO
558
559!
560!--          Calculate tendencies for the next Runge-Kutta step
561             IF ( timestep_scheme(1:5) == 'runge' )  THEN
562                IF ( intermediate_timestep_count == 1 )  THEN
563                   DO  k = nzb+1, nzt
564                      tv_m(k,j,i) = tend(k,j,i)
565                   ENDDO
566                ELSEIF ( intermediate_timestep_count < &
567                         intermediate_timestep_count_max )  THEN
568                   DO  k = nzb+1, nzt
569                      tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
570                                     + 5.3125_wp * tv_m(k,j,i)
571                   ENDDO
572                ENDIF
573             ENDIF
574
575          ENDIF
576
577!
578!--       Tendency terms for w-velocity component
579          tend(:,j,i) = 0.0_wp
580          IF ( timestep_scheme(1:5) == 'runge' )  THEN
581             IF ( ws_scheme_mom )  THEN
582                CALL advec_w_ws( i, j, i_omp_start, tn )
583             ELSE
584                CALL advec_w_pw( i, j )
585             END IF
586          ELSE
587             CALL advec_w_up( i, j )
588          ENDIF
589          CALL diffusion_w( i, j )
590          CALL coriolis( i, j, 3 )
591
592          IF ( .NOT. neutral )  THEN
593             IF ( ocean )  THEN
594                CALL buoyancy( i, j, rho_ocean, 3 )
595             ELSE
596                IF ( .NOT. humidity )  THEN
597                   CALL buoyancy( i, j, pt, 3 )
598                ELSE
599                   CALL buoyancy( i, j, vpt, 3 )
600                ENDIF
601             ENDIF
602          ENDIF
603
604!
605!--       Drag by plant canopy
606          IF ( plant_canopy )  CALL pcm_tendency( i, j, 3 )
607
608!
609!--       Forces by wind turbines
610          IF ( wind_turbine )  CALL wtm_tendencies( i, j, 3 )
611
612          CALL user_actions( i, j, 'w-tendency' )
613!
614!--       Prognostic equation for w-velocity component
615          DO  k = nzb+1, nzt-1
616             w_p(k,j,i) = w(k,j,i) + ( dt_3d *                                 &
617                                             ( tsc(2) * tend(k,j,i) +          &
618                                               tsc(3) * tw_m(k,j,i) )          &
619                                             - tsc(5) * rdf(k) * w(k,j,i)      &
620                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
621                                                BTEST( wall_flags_0(k,j,i), 3 )&
622                                              )
623          ENDDO
624
625!
626!--       Calculate tendencies for the next Runge-Kutta step
627          IF ( timestep_scheme(1:5) == 'runge' )  THEN
628             IF ( intermediate_timestep_count == 1 )  THEN
629                DO  k = nzb+1, nzt-1
630                   tw_m(k,j,i) = tend(k,j,i)
631                ENDDO
632             ELSEIF ( intermediate_timestep_count < &
633                      intermediate_timestep_count_max )  THEN
634                DO  k = nzb+1, nzt-1
635                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
636                                  + 5.3125_wp * tw_m(k,j,i)
637                ENDDO
638             ENDIF
639          ENDIF
640
641!
642!--       If required, compute prognostic equation for potential temperature
643          IF ( .NOT. neutral )  THEN
644!
645!--          Tendency terms for potential temperature
646             tend(:,j,i) = 0.0_wp
647             IF ( timestep_scheme(1:5) == 'runge' )  THEN
648                   IF ( ws_scheme_sca )  THEN
649                      CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, diss_s_pt,   &
650                                       flux_l_pt, diss_l_pt, i_omp_start, tn )
651                   ELSE
652                      CALL advec_s_pw( i, j, pt )
653                   ENDIF
654             ELSE
655                CALL advec_s_up( i, j, pt )
656             ENDIF
657             CALL diffusion_s( i, j, pt,                                       &
658                               surf_def_h(0)%shf, surf_def_h(1)%shf,           &
659                               surf_def_h(2)%shf,                              &
660                               surf_lsm_h%shf,    surf_usm_h%shf,              &
661                               surf_def_v(0)%shf, surf_def_v(1)%shf,           &
662                               surf_def_v(2)%shf, surf_def_v(3)%shf,           &
663                               surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,           &
664                               surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,           &
665                               surf_usm_v(0)%shf, surf_usm_v(1)%shf,           &
666                               surf_usm_v(2)%shf, surf_usm_v(3)%shf )
667!
668!--          If required compute heating/cooling due to long wave radiation
669!--          processes
670             IF ( cloud_top_radiation )  THEN
671                CALL calc_radiation( i, j )
672             ENDIF
673
674!
675!--          Consideration of heat sources within the plant canopy
676             IF ( plant_canopy  .AND.  cthf /= 0.0_wp )  THEN
677                CALL pcm_tendency( i, j, 4 )
678             ENDIF
679
680!
681!--          Large scale advection
682             IF ( large_scale_forcing )  THEN
683                CALL ls_advec( i, j, simulated_time, 'pt' )
684             ENDIF
685
686!
687!--          Nudging
688             IF ( nudging )  CALL nudge( i, j, simulated_time, 'pt' )
689
690!
691!--          If required, compute effect of large-scale subsidence/ascent
692             IF ( large_scale_subsidence  .AND.                                &
693                  .NOT. use_subsidence_tendencies )  THEN
694                CALL subsidence( i, j, tend, pt, pt_init, 2 )
695             ENDIF
696
697!
698!--          If required, add tendency due to radiative heating/cooling
699             IF ( radiation  .AND.                                             &
700                  simulated_time > skip_time_do_radiation )  THEN
701                CALL radiation_tendency ( i, j, tend )
702             ENDIF
703
704
705             CALL user_actions( i, j, 'pt-tendency' )
706!
707!--          Prognostic equation for potential temperature
708             DO  k = nzb+1, nzt
709                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d *                            &
710                                                  ( tsc(2) * tend(k,j,i) +     &
711                                                    tsc(3) * tpt_m(k,j,i) )    &
712                                                  - tsc(5)                     &
713                                                  * ( pt(k,j,i) - pt_init(k) ) &
714                                                  * ( rdf_sc(k) + ptdf_x(i)    &
715                                                                + ptdf_y(j) )  &
716                                          )                                    &
717                                       * MERGE( 1.0_wp, 0.0_wp,                &
718                                                BTEST( wall_flags_0(k,j,i), 0 )&
719                                              )
720             ENDDO
721
722!
723!--          Calculate tendencies for the next Runge-Kutta step
724             IF ( timestep_scheme(1:5) == 'runge' )  THEN
725                IF ( intermediate_timestep_count == 1 )  THEN
726                   DO  k = nzb+1, nzt
727                      tpt_m(k,j,i) = tend(k,j,i)
728                   ENDDO
729                ELSEIF ( intermediate_timestep_count < &
730                         intermediate_timestep_count_max )  THEN
731                   DO  k = nzb+1, nzt
732                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
733                                        5.3125_wp * tpt_m(k,j,i)
734                   ENDDO
735                ENDIF
736             ENDIF
737
738          ENDIF
739
740!
741!--       If required, compute prognostic equation for salinity
742          IF ( ocean )  THEN
743
744!
745!--          Tendency-terms for salinity
746             tend(:,j,i) = 0.0_wp
747             IF ( timestep_scheme(1:5) == 'runge' ) &
748             THEN
749                IF ( ws_scheme_sca )  THEN
750                    CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa,  &
751                                diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn  )
752                ELSE
753                    CALL advec_s_pw( i, j, sa )
754                ENDIF
755             ELSE
756                CALL advec_s_up( i, j, sa )
757             ENDIF
758             CALL diffusion_s( i, j, sa,                                       &
759                               surf_def_h(0)%sasws, surf_def_h(1)%sasws,       &
760                               surf_def_h(2)%sasws,                            &
761                               surf_lsm_h%sasws,    surf_usm_h%sasws,          & 
762                               surf_def_v(0)%sasws, surf_def_v(1)%sasws,       &
763                               surf_def_v(2)%sasws, surf_def_v(3)%sasws,       &
764                               surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws,       &
765                               surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws,       &
766                               surf_usm_v(0)%sasws, surf_usm_v(1)%sasws,       &
767                               surf_usm_v(2)%sasws, surf_usm_v(3)%sasws )
768
769             CALL user_actions( i, j, 'sa-tendency' )
770
771!
772!--          Prognostic equation for salinity
773             DO  k = nzb+1, nzt
774                sa_p(k,j,i) = sa(k,j,i) + ( dt_3d *                            &
775                                                  ( tsc(2) * tend(k,j,i) +     &
776                                                    tsc(3) * tsa_m(k,j,i) )    &
777                                                  - tsc(5) * rdf_sc(k)         &
778                                                   * ( sa(k,j,i) - sa_init(k) )&
779                                          ) * MERGE( 1.0_wp, 0.0_wp,           &
780                                                BTEST( wall_flags_0(k,j,i), 0 )&
781                                                   )
782                IF ( sa_p(k,j,i) < 0.0_wp )  sa_p(k,j,i) = 0.1_wp * sa(k,j,i)
783             ENDDO
784
785!
786!--          Calculate tendencies for the next Runge-Kutta step
787             IF ( timestep_scheme(1:5) == 'runge' )  THEN
788                IF ( intermediate_timestep_count == 1 )  THEN
789                   DO  k = nzb+1, nzt
790                      tsa_m(k,j,i) = tend(k,j,i)
791                   ENDDO
792                ELSEIF ( intermediate_timestep_count < &
793                         intermediate_timestep_count_max )  THEN
794                   DO  k = nzb+1, nzt
795                      tsa_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
796                                        5.3125_wp * tsa_m(k,j,i)
797                   ENDDO 
798                ENDIF
799             ENDIF
800
801!
802!--          Calculate density by the equation of state for seawater
803             CALL eqn_state_seawater( i, j )
804
805          ENDIF
806
807!
808!--       If required, compute prognostic equation for total water content
809          IF ( humidity )  THEN
810
811!
812!--          Tendency-terms for total water content / scalar
813             tend(:,j,i) = 0.0_wp
814             IF ( timestep_scheme(1:5) == 'runge' ) &
815             THEN
816                IF ( ws_scheme_sca )  THEN
817                   CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
818                                diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
819                ELSE
820                   CALL advec_s_pw( i, j, q )
821                ENDIF
822             ELSE
823                CALL advec_s_up( i, j, q )
824             ENDIF
825             CALL diffusion_s( i, j, q,                                        &
826                               surf_def_h(0)%qsws, surf_def_h(1)%qsws,         &
827                               surf_def_h(2)%qsws,                             &
828                               surf_lsm_h%qsws,    surf_usm_h%qsws,            &
829                               surf_def_v(0)%qsws, surf_def_v(1)%qsws,         &
830                               surf_def_v(2)%qsws, surf_def_v(3)%qsws,         &
831                               surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,         &
832                               surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,         &
833                               surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,         &
834                               surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
835
836!
837!--          Sink or source of humidity due to canopy elements
838             IF ( plant_canopy )  CALL pcm_tendency( i, j, 5 )
839
840!
841!--          Large scale advection
842             IF ( large_scale_forcing )  THEN
843                CALL ls_advec( i, j, simulated_time, 'q' )
844             ENDIF
845
846!
847!--          Nudging
848             IF ( nudging )  CALL nudge( i, j, simulated_time, 'q' )
849
850!
851!--          If required compute influence of large-scale subsidence/ascent
852             IF ( large_scale_subsidence  .AND.                                &
853                  .NOT. use_subsidence_tendencies )  THEN
854                CALL subsidence( i, j, tend, q, q_init, 3 )
855             ENDIF
856
857             CALL user_actions( i, j, 'q-tendency' )
858
859!
860!--          Prognostic equation for total water content / scalar
861             DO  k = nzb+1, nzt
862                q_p(k,j,i) = q(k,j,i) + ( dt_3d *                              &
863                                                ( tsc(2) * tend(k,j,i) +       &
864                                                  tsc(3) * tq_m(k,j,i) )       &
865                                                - tsc(5) * rdf_sc(k) *         &
866                                                      ( q(k,j,i) - q_init(k) ) &
867                                        )                                      &
868                                       * MERGE( 1.0_wp, 0.0_wp,                &
869                                                BTEST( wall_flags_0(k,j,i), 0 )&
870                                              )               
871                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
872             ENDDO
873
874!
875!--          Calculate tendencies for the next Runge-Kutta step
876             IF ( timestep_scheme(1:5) == 'runge' )  THEN
877                IF ( intermediate_timestep_count == 1 )  THEN
878                   DO  k = nzb+1, nzt
879                      tq_m(k,j,i) = tend(k,j,i)
880                   ENDDO
881                ELSEIF ( intermediate_timestep_count < &
882                         intermediate_timestep_count_max )  THEN
883                   DO  k = nzb+1, nzt
884                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
885                                       5.3125_wp * tq_m(k,j,i)
886                   ENDDO
887                ENDIF
888             ENDIF
889
890!
891!--          If required, calculate prognostic equations for rain water content
892!--          and rain drop concentration
893             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
894!
895!--             Calculate prognostic equation for rain water content
896                tend(:,j,i) = 0.0_wp
897                IF ( timestep_scheme(1:5) == 'runge' ) &
898                THEN
899                   IF ( ws_scheme_sca )  THEN
900                      CALL advec_s_ws( i, j, qr, 'qr', flux_s_qr,       &
901                                       diss_s_qr, flux_l_qr, diss_l_qr, &
902                                       i_omp_start, tn )
903                   ELSE
904                      CALL advec_s_pw( i, j, qr )
905                   ENDIF
906                ELSE
907                   CALL advec_s_up( i, j, qr )
908                ENDIF
909                CALL diffusion_s( i, j, qr,                                   &
910                                  surf_def_h(0)%qrsws, surf_def_h(1)%qrsws,   &
911                                  surf_def_h(2)%qrsws,                        &
912                                  surf_lsm_h%qrsws,    surf_usm_h%qrsws,      & 
913                                  surf_def_v(0)%qrsws, surf_def_v(1)%qrsws,   &
914                                  surf_def_v(2)%qrsws, surf_def_v(3)%qrsws,   &
915                                  surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws,   &
916                                  surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws,   &
917                                  surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws,   &
918                                  surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws )
919
920!
921!--             Prognostic equation for rain water content
922                DO  k = nzb+1, nzt
923                   qr_p(k,j,i) = qr(k,j,i) + ( dt_3d *                         &
924                                                      ( tsc(2) * tend(k,j,i) + &
925                                                        tsc(3) * tqr_m(k,j,i) )&
926                                                      - tsc(5) * rdf_sc(k)     &
927                                                               * qr(k,j,i)     &
928                                             )                                 &
929                                       * MERGE( 1.0_wp, 0.0_wp,                &
930                                                BTEST( wall_flags_0(k,j,i), 0 )&
931                                              ) 
932                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
933                ENDDO
934!
935!--             Calculate tendencies for the next Runge-Kutta step
936                IF ( timestep_scheme(1:5) == 'runge' )  THEN
937                   IF ( intermediate_timestep_count == 1 )  THEN
938                      DO  k = nzb+1, nzt
939                         tqr_m(k,j,i) = tend(k,j,i)
940                      ENDDO
941                   ELSEIF ( intermediate_timestep_count < &
942                            intermediate_timestep_count_max )  THEN
943                      DO  k = nzb+1, nzt
944                         tqr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
945                                           5.3125_wp * tqr_m(k,j,i)
946                      ENDDO
947                   ENDIF
948                ENDIF
949
950!
951!--             Calculate prognostic equation for rain drop concentration.
952                tend(:,j,i) = 0.0_wp
953                IF ( timestep_scheme(1:5) == 'runge' )  THEN
954                   IF ( ws_scheme_sca )  THEN
955                      CALL advec_s_ws( i, j, nr, 'nr', flux_s_nr,    &
956                                    diss_s_nr, flux_l_nr, diss_l_nr, &
957                                    i_omp_start, tn )
958                   ELSE
959                      CALL advec_s_pw( i, j, nr )
960                   ENDIF
961                ELSE
962                   CALL advec_s_up( i, j, nr )
963                ENDIF
964                CALL diffusion_s( i, j, nr,                                    &
965                                  surf_def_h(0)%nrsws, surf_def_h(1)%nrsws,    &
966                                  surf_def_h(2)%nrsws,                         &
967                                  surf_lsm_h%nrsws,    surf_usm_h%nrsws,       &
968                                  surf_def_v(0)%nrsws, surf_def_v(1)%nrsws,    &
969                                  surf_def_v(2)%nrsws, surf_def_v(3)%nrsws,    &
970                                  surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws,    &
971                                  surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws,    &
972                                  surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws,    &
973                                  surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws )
974
975!
976!--             Prognostic equation for rain drop concentration
977                DO  k = nzb+1, nzt
978                   nr_p(k,j,i) = nr(k,j,i) + ( dt_3d *                         &
979                                                      ( tsc(2) * tend(k,j,i) + &
980                                                        tsc(3) * tnr_m(k,j,i) )&
981                                                      - tsc(5) * rdf_sc(k)     &
982                                                               * nr(k,j,i)     &
983                                             )                                 &
984                                       * MERGE( 1.0_wp, 0.0_wp,                &
985                                                BTEST( wall_flags_0(k,j,i), 0 )&
986                                              )
987                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
988                ENDDO
989!
990!--             Calculate tendencies for the next Runge-Kutta step
991                IF ( timestep_scheme(1:5) == 'runge' )  THEN
992                   IF ( intermediate_timestep_count == 1 )  THEN
993                      DO  k = nzb+1, nzt
994                         tnr_m(k,j,i) = tend(k,j,i)
995                      ENDDO
996                   ELSEIF ( intermediate_timestep_count < &
997                            intermediate_timestep_count_max )  THEN
998                      DO  k = nzb+1, nzt
999                         tnr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
1000                                           5.3125_wp * tnr_m(k,j,i)
1001                      ENDDO
1002                   ENDIF
1003                ENDIF
1004
1005             ENDIF
1006
1007          ENDIF
1008
1009!
1010!--       If required, compute prognostic equation for scalar
1011          IF ( passive_scalar )  THEN
1012!
1013!--          Tendency-terms for total water content / scalar
1014             tend(:,j,i) = 0.0_wp
1015             IF ( timestep_scheme(1:5) == 'runge' ) &
1016             THEN
1017                IF ( ws_scheme_sca )  THEN
1018                   CALL advec_s_ws( i, j, s, 's', flux_s_s, &
1019                                diss_s_s, flux_l_s, diss_l_s, i_omp_start, tn )
1020                ELSE
1021                   CALL advec_s_pw( i, j, s )
1022                ENDIF
1023             ELSE
1024                CALL advec_s_up( i, j, s )
1025             ENDIF
1026             CALL diffusion_s( i, j, s,                                        &
1027                               surf_def_h(0)%ssws, surf_def_h(1)%ssws,         &
1028                               surf_def_h(2)%ssws,                             &
1029                               surf_lsm_h%ssws,    surf_usm_h%ssws,            &
1030                               surf_def_v(0)%ssws, surf_def_v(1)%ssws,         &
1031                               surf_def_v(2)%ssws, surf_def_v(3)%ssws,         &
1032                               surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,         &
1033                               surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,         &
1034                               surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,         &
1035                               surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1036
1037!
1038!--          Sink or source of scalar concentration due to canopy elements
1039             IF ( plant_canopy )  CALL pcm_tendency( i, j, 7 )
1040
1041!
1042!--          Large scale advection, still need to be extended for scalars
1043!              IF ( large_scale_forcing )  THEN
1044!                 CALL ls_advec( i, j, simulated_time, 's' )
1045!              ENDIF
1046
1047!
1048!--          Nudging, still need to be extended for scalars
1049!              IF ( nudging )  CALL nudge( i, j, simulated_time, 's' )
1050
1051!
1052!--          If required compute influence of large-scale subsidence/ascent.
1053!--          Note, the last argument is of no meaning in this case, as it is
1054!--          only used in conjunction with large_scale_forcing, which is to
1055!--          date not implemented for scalars.
1056             IF ( large_scale_subsidence  .AND.                                &
1057                  .NOT. use_subsidence_tendencies  .AND.                       &
1058                  .NOT. large_scale_forcing )  THEN
1059                CALL subsidence( i, j, tend, s, s_init, 3 )
1060             ENDIF
1061
1062             CALL user_actions( i, j, 's-tendency' )
1063
1064!
1065!--          Prognostic equation for scalar
1066             DO  k = nzb+1, nzt
1067                s_p(k,j,i) = s(k,j,i) + (  dt_3d *                             &
1068                                            ( tsc(2) * tend(k,j,i) +           &
1069                                              tsc(3) * ts_m(k,j,i) )           &
1070                                            - tsc(5) * rdf_sc(k)               &
1071                                                     * ( s(k,j,i) - s_init(k) )&
1072                                        )                                      &
1073                                       * MERGE( 1.0_wp, 0.0_wp,                &
1074                                                BTEST( wall_flags_0(k,j,i), 0 )&
1075                                              )
1076                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
1077             ENDDO
1078
1079!
1080!--          Calculate tendencies for the next Runge-Kutta step
1081             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1082                IF ( intermediate_timestep_count == 1 )  THEN
1083                   DO  k = nzb+1, nzt
1084                      ts_m(k,j,i) = tend(k,j,i)
1085                   ENDDO
1086                ELSEIF ( intermediate_timestep_count < &
1087                         intermediate_timestep_count_max )  THEN
1088                   DO  k = nzb+1, nzt
1089                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
1090                                       5.3125_wp * ts_m(k,j,i)
1091                   ENDDO
1092                ENDIF
1093             ENDIF
1094
1095          ENDIF
1096!
1097!--       If required, compute prognostic equation for turbulent kinetic
1098!--       energy (TKE)
1099          IF ( .NOT. constant_diffusion )  THEN
1100
1101!
1102!--          Tendency-terms for TKE
1103             tend(:,j,i) = 0.0_wp
1104             IF ( timestep_scheme(1:5) == 'runge'  &
1105                 .AND.  .NOT. use_upstream_for_tke )  THEN
1106                 IF ( ws_scheme_sca )  THEN
1107                     CALL advec_s_ws( i, j, e, 'e', flux_s_e, diss_s_e, &
1108                                      flux_l_e, diss_l_e , i_omp_start, tn )
1109                 ELSE
1110                     CALL advec_s_pw( i, j, e )
1111                 ENDIF
1112             ELSE
1113                CALL advec_s_up( i, j, e )
1114             ENDIF
1115             IF ( .NOT. humidity )  THEN
1116                IF ( ocean )  THEN
1117                   CALL diffusion_e( i, j, prho, prho_reference )
1118                ELSE
1119                   CALL diffusion_e( i, j, pt, pt_reference )
1120                ENDIF
1121             ELSE
1122                CALL diffusion_e( i, j, vpt, pt_reference )
1123             ENDIF
1124             CALL production_e( i, j )
1125
1126!
1127!--          Additional sink term for flows through plant canopies
1128             IF ( plant_canopy )  CALL pcm_tendency( i, j, 6 )
1129
1130             CALL user_actions( i, j, 'e-tendency' )
1131
1132!
1133!--          Prognostic equation for TKE.
1134!--          Eliminate negative TKE values, which can occur due to numerical
1135!--          reasons in the course of the integration. In such cases the old
1136!--          TKE value is reduced by 90%.
1137             DO  k = nzb+1, nzt
1138                e_p(k,j,i) = e(k,j,i) + ( dt_3d *                              &
1139                                                ( tsc(2) * tend(k,j,i) +       &
1140                                                  tsc(3) * te_m(k,j,i) )       &
1141                                        )                                      &
1142                                       * MERGE( 1.0_wp, 0.0_wp,                &
1143                                                BTEST( wall_flags_0(k,j,i), 0 )&
1144                                              )
1145                IF ( e_p(k,j,i) < 0.0_wp )  e_p(k,j,i) = 0.1_wp * e(k,j,i)
1146             ENDDO
1147
1148!
1149!--          Calculate tendencies for the next Runge-Kutta step
1150             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1151                IF ( intermediate_timestep_count == 1 )  THEN
1152                   DO  k = nzb+1, nzt
1153                      te_m(k,j,i) = tend(k,j,i)
1154                   ENDDO
1155                ELSEIF ( intermediate_timestep_count < &
1156                         intermediate_timestep_count_max )  THEN
1157                   DO  k = nzb+1, nzt
1158                      te_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
1159                                       5.3125_wp * te_m(k,j,i)
1160                   ENDDO
1161                ENDIF
1162             ENDIF
1163
1164          ENDIF   ! TKE equation
1165
1166       ENDDO
1167    ENDDO
1168    !$OMP END PARALLEL
1169
1170
1171
1172
1173    CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1174
1175
1176 END SUBROUTINE prognostic_equations_cache
1177
1178
1179!------------------------------------------------------------------------------!
1180! Description:
1181! ------------
1182!> Version for vector machines
1183!------------------------------------------------------------------------------!
1184
1185 SUBROUTINE prognostic_equations_vector
1186
1187
1188    IMPLICIT NONE
1189
1190    INTEGER(iwp) ::  i    !<
1191    INTEGER(iwp) ::  j    !<
1192    INTEGER(iwp) ::  k    !<
1193
1194    REAL(wp)     ::  sbt  !<
1195
1196
1197!
1198!-- If required, calculate cloud microphysical impacts
1199    IF ( cloud_physics  .AND.  .NOT. microphysics_sat_adjust  .AND.            &
1200         ( intermediate_timestep_count == 1  .OR.                              &
1201           call_microphysics_at_all_substeps )                                 &
1202       )  THEN
1203       CALL cpu_log( log_point(51), 'microphysics', 'start' )
1204       CALL microphysics_control
1205       CALL cpu_log( log_point(51), 'microphysics', 'stop' )
1206    ENDIF
1207
1208!
1209!-- u-velocity component
1210    CALL cpu_log( log_point(5), 'u-equation', 'start' )
1211
1212    tend = 0.0_wp
1213    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1214       IF ( ws_scheme_mom )  THEN
1215          CALL advec_u_ws
1216       ELSE
1217          CALL advec_u_pw
1218       ENDIF
1219    ELSE
1220       CALL advec_u_up
1221    ENDIF
1222    CALL diffusion_u
1223    CALL coriolis( 1 )
1224    IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
1225       CALL buoyancy( pt, 1 )
1226    ENDIF
1227
1228!
1229!-- Drag by plant canopy
1230    IF ( plant_canopy )  CALL pcm_tendency( 1 )
1231
1232!
1233!-- External pressure gradient
1234    IF ( dp_external )  THEN
1235       DO  i = nxlu, nxr
1236          DO  j = nys, nyn
1237             DO  k = dp_level_ind_b+1, nzt
1238                tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
1239             ENDDO
1240          ENDDO
1241       ENDDO
1242    ENDIF
1243
1244!
1245!-- Nudging
1246    IF ( nudging )  CALL nudge( simulated_time, 'u' )
1247
1248!
1249!-- Forces by wind turbines
1250    IF ( wind_turbine )  CALL wtm_tendencies( 1 )
1251
1252    CALL user_actions( 'u-tendency' )
1253
1254!
1255!-- Prognostic equation for u-velocity component
1256    DO  i = nxlu, nxr
1257       DO  j = nys, nyn
1258          DO  k = nzb+1, nzt
1259             u_p(k,j,i) = u(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +          &
1260                                                 tsc(3) * tu_m(k,j,i) )          &
1261                                               - tsc(5) * rdf(k) *               &
1262                                                        ( u(k,j,i) - u_init(k) ) &
1263                                     ) * MERGE( 1.0_wp, 0.0_wp,                  &
1264                                                 BTEST( wall_flags_0(k,j,i), 1 ) &
1265                                              )
1266          ENDDO
1267       ENDDO
1268    ENDDO
1269
1270!
1271!-- Calculate tendencies for the next Runge-Kutta step
1272    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1273       IF ( intermediate_timestep_count == 1 )  THEN
1274          DO  i = nxlu, nxr
1275             DO  j = nys, nyn
1276                DO  k = nzb+1, nzt
1277                   tu_m(k,j,i) = tend(k,j,i)
1278                ENDDO
1279             ENDDO
1280          ENDDO
1281       ELSEIF ( intermediate_timestep_count < &
1282                intermediate_timestep_count_max )  THEN
1283          DO  i = nxlu, nxr
1284             DO  j = nys, nyn
1285                DO  k = nzb+1, nzt
1286                   tu_m(k,j,i) =    -9.5625_wp * tend(k,j,i)                   &
1287                                   + 5.3125_wp * tu_m(k,j,i)
1288                ENDDO
1289             ENDDO
1290          ENDDO
1291       ENDIF
1292    ENDIF
1293
1294    CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1295
1296!
1297!-- v-velocity component
1298    CALL cpu_log( log_point(6), 'v-equation', 'start' )
1299
1300    tend = 0.0_wp
1301    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1302       IF ( ws_scheme_mom )  THEN
1303          CALL advec_v_ws
1304       ELSE
1305          CALL advec_v_pw
1306       END IF
1307    ELSE
1308       CALL advec_v_up
1309    ENDIF
1310    CALL diffusion_v
1311    CALL coriolis( 2 )
1312
1313!
1314!-- Drag by plant canopy
1315    IF ( plant_canopy )  CALL pcm_tendency( 2 )
1316
1317!
1318!-- External pressure gradient
1319    IF ( dp_external )  THEN
1320       DO  i = nxl, nxr
1321          DO  j = nysv, nyn
1322             DO  k = dp_level_ind_b+1, nzt
1323                tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1324             ENDDO
1325          ENDDO
1326       ENDDO
1327    ENDIF
1328
1329!
1330!-- Nudging
1331    IF ( nudging )  CALL nudge( simulated_time, 'v' )
1332
1333!
1334!-- Forces by wind turbines
1335    IF ( wind_turbine )  CALL wtm_tendencies( 2 )
1336
1337    CALL user_actions( 'v-tendency' )
1338
1339!
1340!-- Prognostic equation for v-velocity component
1341    DO  i = nxl, nxr
1342       DO  j = nysv, nyn
1343          DO  k = nzb+1, nzt
1344             v_p(k,j,i) = v(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1345                                                 tsc(3) * tv_m(k,j,i) )        &
1346                                               - tsc(5) * rdf(k) *             &
1347                                                      ( v(k,j,i) - v_init(k) ) &
1348                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1349                                                BTEST( wall_flags_0(k,j,i), 2 )&
1350                                              )
1351          ENDDO
1352       ENDDO
1353    ENDDO
1354
1355!
1356!-- Calculate tendencies for the next Runge-Kutta step
1357    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1358       IF ( intermediate_timestep_count == 1 )  THEN
1359          DO  i = nxl, nxr
1360             DO  j = nysv, nyn
1361                DO  k = nzb+1, nzt
1362                   tv_m(k,j,i) = tend(k,j,i)
1363                ENDDO
1364             ENDDO
1365          ENDDO
1366       ELSEIF ( intermediate_timestep_count < &
1367                intermediate_timestep_count_max )  THEN
1368          DO  i = nxl, nxr
1369             DO  j = nysv, nyn
1370                DO  k = nzb+1, nzt
1371                   tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1372                                  + 5.3125_wp * tv_m(k,j,i)
1373                ENDDO
1374             ENDDO
1375          ENDDO
1376       ENDIF
1377    ENDIF
1378
1379    CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1380
1381!
1382!-- w-velocity component
1383    CALL cpu_log( log_point(7), 'w-equation', 'start' )
1384
1385    tend = 0.0_wp
1386    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1387       IF ( ws_scheme_mom )  THEN
1388          CALL advec_w_ws
1389       ELSE
1390          CALL advec_w_pw
1391       ENDIF
1392    ELSE
1393       CALL advec_w_up
1394    ENDIF
1395    CALL diffusion_w
1396    CALL coriolis( 3 )
1397
1398    IF ( .NOT. neutral )  THEN
1399       IF ( ocean )  THEN
1400          CALL buoyancy( rho_ocean, 3 )
1401       ELSE
1402          IF ( .NOT. humidity )  THEN
1403             CALL buoyancy( pt, 3 )
1404          ELSE
1405             CALL buoyancy( vpt, 3 )
1406          ENDIF
1407       ENDIF
1408    ENDIF
1409
1410!
1411!-- Drag by plant canopy
1412    IF ( plant_canopy )  CALL pcm_tendency( 3 )
1413
1414!
1415!-- Forces by wind turbines
1416    IF ( wind_turbine )  CALL wtm_tendencies( 3 )
1417
1418    CALL user_actions( 'w-tendency' )
1419
1420!
1421!-- Prognostic equation for w-velocity component
1422    DO  i = nxl, nxr
1423       DO  j = nys, nyn
1424          DO  k = nzb+1, nzt-1
1425             w_p(k,j,i) = w(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1426                                                 tsc(3) * tw_m(k,j,i) )        &
1427                                               - tsc(5) * rdf(k) * w(k,j,i)    &
1428                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1429                                                BTEST( wall_flags_0(k,j,i), 3 )&
1430                                              )
1431          ENDDO
1432       ENDDO
1433    ENDDO
1434
1435!
1436!-- Calculate tendencies for the next Runge-Kutta step
1437    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1438       IF ( intermediate_timestep_count == 1 )  THEN
1439          DO  i = nxl, nxr
1440             DO  j = nys, nyn
1441                DO  k = nzb+1, nzt-1
1442                   tw_m(k,j,i) = tend(k,j,i)
1443                ENDDO
1444             ENDDO
1445          ENDDO
1446       ELSEIF ( intermediate_timestep_count < &
1447                intermediate_timestep_count_max )  THEN
1448          DO  i = nxl, nxr
1449             DO  j = nys, nyn
1450                DO  k = nzb+1, nzt-1
1451                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1452                                  + 5.3125_wp * tw_m(k,j,i)
1453                ENDDO
1454             ENDDO
1455          ENDDO
1456       ENDIF
1457    ENDIF
1458
1459    CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1460
1461
1462!
1463!-- If required, compute prognostic equation for potential temperature
1464    IF ( .NOT. neutral )  THEN
1465
1466       CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1467
1468!
1469!--    pt-tendency terms with communication
1470       sbt = tsc(2)
1471       IF ( scalar_advec == 'bc-scheme' )  THEN
1472
1473          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1474!
1475!--          Bott-Chlond scheme always uses Euler time step. Thus:
1476             sbt = 1.0_wp
1477          ENDIF
1478          tend = 0.0_wp
1479          CALL advec_s_bc( pt, 'pt' )
1480
1481       ENDIF
1482
1483!
1484!--    pt-tendency terms with no communication
1485       IF ( scalar_advec /= 'bc-scheme' )  THEN
1486          tend = 0.0_wp
1487          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1488             IF ( ws_scheme_sca )  THEN
1489                CALL advec_s_ws( pt, 'pt' )
1490             ELSE
1491                CALL advec_s_pw( pt )
1492             ENDIF
1493          ELSE
1494             CALL advec_s_up( pt )
1495          ENDIF
1496       ENDIF
1497
1498       CALL diffusion_s( pt,                                                   &
1499                         surf_def_h(0)%shf, surf_def_h(1)%shf,                 &
1500                         surf_def_h(2)%shf,                                    &
1501                         surf_lsm_h%shf,    surf_usm_h%shf,                    &
1502                         surf_def_v(0)%shf, surf_def_v(1)%shf,                 &
1503                         surf_def_v(2)%shf, surf_def_v(3)%shf,                 &
1504                         surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,                 &
1505                         surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,                 &
1506                         surf_usm_v(0)%shf, surf_usm_v(1)%shf,                 &
1507                         surf_usm_v(2)%shf, surf_usm_v(3)%shf )
1508!
1509!--    If required compute heating/cooling due to long wave radiation processes
1510       IF ( cloud_top_radiation )  THEN
1511          CALL calc_radiation
1512       ENDIF
1513
1514!
1515!--    Consideration of heat sources within the plant canopy
1516       IF ( plant_canopy .AND. ( cthf /= 0.0_wp ) ) THEN
1517          CALL pcm_tendency( 4 )
1518       ENDIF
1519
1520!
1521!--    Large scale advection
1522       IF ( large_scale_forcing )  THEN
1523          CALL ls_advec( simulated_time, 'pt' )
1524       ENDIF
1525
1526!
1527!--    Nudging
1528       IF ( nudging )  CALL nudge( simulated_time, 'pt' )
1529
1530!
1531!--    If required compute influence of large-scale subsidence/ascent
1532       IF ( large_scale_subsidence  .AND.                                      &
1533            .NOT. use_subsidence_tendencies )  THEN
1534          CALL subsidence( tend, pt, pt_init, 2 )
1535       ENDIF
1536
1537!
1538!--    If required, add tendency due to radiative heating/cooling
1539       IF ( radiation  .AND.                                                   &
1540            simulated_time > skip_time_do_radiation )  THEN
1541            CALL radiation_tendency ( tend )
1542       ENDIF
1543
1544       CALL user_actions( 'pt-tendency' )
1545
1546!
1547!--    Prognostic equation for potential temperature
1548       DO  i = nxl, nxr
1549          DO  j = nys, nyn
1550             DO  k = nzb+1, nzt
1551                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +      &
1552                                                   tsc(3) * tpt_m(k,j,i) )     &
1553                                                 - tsc(5) *                    &
1554                                                   ( pt(k,j,i) - pt_init(k) ) *&
1555                                          ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )&
1556                                          )                                    &
1557                                   * MERGE( 1.0_wp, 0.0_wp,                    &
1558                                             BTEST( wall_flags_0(k,j,i), 0 )   &
1559                                          )
1560             ENDDO
1561          ENDDO
1562       ENDDO
1563
1564!
1565!--    Calculate tendencies for the next Runge-Kutta step
1566       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1567          IF ( intermediate_timestep_count == 1 )  THEN
1568             DO  i = nxl, nxr
1569                DO  j = nys, nyn
1570                   DO  k = nzb+1, nzt
1571                      tpt_m(k,j,i) = tend(k,j,i)
1572                   ENDDO
1573                ENDDO
1574             ENDDO
1575          ELSEIF ( intermediate_timestep_count < &
1576                   intermediate_timestep_count_max )  THEN
1577             DO  i = nxl, nxr
1578                DO  j = nys, nyn
1579                   DO  k = nzb+1, nzt
1580                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
1581                                        5.3125_wp * tpt_m(k,j,i)
1582                   ENDDO
1583                ENDDO
1584             ENDDO
1585          ENDIF
1586       ENDIF
1587
1588       CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1589
1590    ENDIF
1591
1592!
1593!-- If required, compute prognostic equation for salinity
1594    IF ( ocean )  THEN
1595
1596       CALL cpu_log( log_point(37), 'sa-equation', 'start' )
1597
1598!
1599!--    sa-tendency terms with communication
1600       sbt = tsc(2)
1601       IF ( scalar_advec == 'bc-scheme' )  THEN
1602
1603          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1604!
1605!--          Bott-Chlond scheme always uses Euler time step. Thus:
1606             sbt = 1.0_wp
1607          ENDIF
1608          tend = 0.0_wp
1609          CALL advec_s_bc( sa, 'sa' )
1610
1611       ENDIF
1612
1613!
1614!--    sa-tendency terms with no communication
1615       IF ( scalar_advec /= 'bc-scheme' )  THEN
1616          tend = 0.0_wp
1617          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1618             IF ( ws_scheme_sca )  THEN
1619                 CALL advec_s_ws( sa, 'sa' )
1620             ELSE
1621                 CALL advec_s_pw( sa )
1622             ENDIF
1623          ELSE
1624             CALL advec_s_up( sa )
1625          ENDIF
1626       ENDIF
1627
1628       CALL diffusion_s( sa,                                                   &
1629                         surf_def_h(0)%sasws, surf_def_h(1)%sasws,             &
1630                         surf_def_h(2)%sasws,                                  &
1631                         surf_lsm_h%sasws,    surf_usm_h%sasws,                &
1632                         surf_def_v(0)%sasws, surf_def_v(1)%sasws,             &
1633                         surf_def_v(2)%sasws, surf_def_v(3)%sasws,             &
1634                         surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws,             &
1635                         surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws,             &
1636                         surf_usm_v(0)%sasws, surf_usm_v(1)%sasws,             &
1637                         surf_usm_v(2)%sasws, surf_usm_v(3)%sasws )
1638
1639       CALL user_actions( 'sa-tendency' )
1640
1641!
1642!--    Prognostic equation for salinity
1643       DO  i = nxl, nxr
1644          DO  j = nys, nyn
1645             DO  k = nzb+1, nzt
1646                sa_p(k,j,i) = sa(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +      &
1647                                                   tsc(3) * tsa_m(k,j,i) )     &
1648                                                 - tsc(5) * rdf_sc(k) *        &
1649                                                 ( sa(k,j,i) - sa_init(k) )    &
1650                                          )                                    &
1651                                   * MERGE( 1.0_wp, 0.0_wp,                    &
1652                                             BTEST( wall_flags_0(k,j,i), 0 )   &
1653                                          )
1654                IF ( sa_p(k,j,i) < 0.0_wp )  sa_p(k,j,i) = 0.1_wp * sa(k,j,i)
1655             ENDDO
1656          ENDDO
1657       ENDDO
1658
1659!
1660!--    Calculate tendencies for the next Runge-Kutta step
1661       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1662          IF ( intermediate_timestep_count == 1 )  THEN
1663             DO  i = nxl, nxr
1664                DO  j = nys, nyn
1665                   DO  k = nzb+1, nzt
1666                      tsa_m(k,j,i) = tend(k,j,i)
1667                   ENDDO
1668                ENDDO
1669             ENDDO
1670          ELSEIF ( intermediate_timestep_count < &
1671                   intermediate_timestep_count_max )  THEN
1672             DO  i = nxl, nxr
1673                DO  j = nys, nyn
1674                   DO  k = nzb+1, nzt
1675                      tsa_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
1676                                        5.3125_wp * tsa_m(k,j,i)
1677                   ENDDO
1678                ENDDO
1679             ENDDO
1680          ENDIF
1681       ENDIF
1682
1683       CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
1684
1685!
1686!--    Calculate density by the equation of state for seawater
1687       CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
1688       CALL eqn_state_seawater
1689       CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
1690
1691    ENDIF
1692
1693!
1694!-- If required, compute prognostic equation for total water content
1695    IF ( humidity )  THEN
1696
1697       CALL cpu_log( log_point(29), 'q-equation', 'start' )
1698
1699!
1700!--    Scalar/q-tendency terms with communication
1701       sbt = tsc(2)
1702       IF ( scalar_advec == 'bc-scheme' )  THEN
1703
1704          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1705!
1706!--          Bott-Chlond scheme always uses Euler time step. Thus:
1707             sbt = 1.0_wp
1708          ENDIF
1709          tend = 0.0_wp
1710          CALL advec_s_bc( q, 'q' )
1711
1712       ENDIF
1713
1714!
1715!--    Scalar/q-tendency terms with no communication
1716       IF ( scalar_advec /= 'bc-scheme' )  THEN
1717          tend = 0.0_wp
1718          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1719             IF ( ws_scheme_sca )  THEN
1720                CALL advec_s_ws( q, 'q' )
1721             ELSE
1722                CALL advec_s_pw( q )
1723             ENDIF
1724          ELSE
1725             CALL advec_s_up( q )
1726          ENDIF
1727       ENDIF
1728
1729       CALL diffusion_s( q,                                                    &
1730                         surf_def_h(0)%qsws, surf_def_h(1)%qsws,               &
1731                         surf_def_h(2)%qsws,                                   &
1732                         surf_lsm_h%qsws,    surf_usm_h%qsws,                  &
1733                         surf_def_v(0)%qsws, surf_def_v(1)%qsws,               &
1734                         surf_def_v(2)%qsws, surf_def_v(3)%qsws,               &
1735                         surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,               &
1736                         surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,               &
1737                         surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,               &
1738                         surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
1739
1740!
1741!--    Sink or source of humidity due to canopy elements
1742       IF ( plant_canopy ) CALL pcm_tendency( 5 )
1743
1744!
1745!--    Large scale advection
1746       IF ( large_scale_forcing )  THEN
1747          CALL ls_advec( simulated_time, 'q' )
1748       ENDIF
1749
1750!
1751!--    Nudging
1752       IF ( nudging )  CALL nudge( simulated_time, 'q' )
1753
1754!
1755!--    If required compute influence of large-scale subsidence/ascent
1756       IF ( large_scale_subsidence  .AND.                                      &
1757            .NOT. use_subsidence_tendencies )  THEN
1758         CALL subsidence( tend, q, q_init, 3 )
1759       ENDIF
1760
1761       CALL user_actions( 'q-tendency' )
1762
1763!
1764!--    Prognostic equation for total water content
1765       DO  i = nxl, nxr
1766          DO  j = nys, nyn
1767             DO  k = nzb+1, nzt
1768                q_p(k,j,i) = q(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
1769                                                 tsc(3) * tq_m(k,j,i) )        &
1770                                               - tsc(5) * rdf_sc(k) *          &
1771                                                      ( q(k,j,i) - q_init(k) ) &
1772                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
1773                                                   BTEST( wall_flags_0(k,j,i), 0 ) &
1774                                                 )
1775                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
1776             ENDDO
1777          ENDDO
1778       ENDDO
1779
1780!
1781!--    Calculate tendencies for the next Runge-Kutta step
1782       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1783          IF ( intermediate_timestep_count == 1 )  THEN
1784             DO  i = nxl, nxr
1785                DO  j = nys, nyn
1786                   DO  k = nzb+1, nzt
1787                      tq_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             DO  i = nxl, nxr
1794                DO  j = nys, nyn
1795                   DO  k = nzb+1, nzt
1796                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
1797                                     + 5.3125_wp * tq_m(k,j,i)
1798                   ENDDO
1799                ENDDO
1800             ENDDO
1801          ENDIF
1802       ENDIF
1803
1804       CALL cpu_log( log_point(29), 'q-equation', 'stop' )
1805
1806!
1807!--    If required, calculate prognostic equations for rain water content
1808!--    and rain drop concentration
1809       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
1810
1811          CALL cpu_log( log_point(52), 'qr-equation', 'start' )
1812
1813!
1814!--       Calculate prognostic equation for rain water content
1815          sbt = tsc(2)
1816          IF ( scalar_advec == 'bc-scheme' )  THEN
1817
1818             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1819!
1820!--             Bott-Chlond scheme always uses Euler time step. Thus:
1821                sbt = 1.0_wp
1822             ENDIF
1823             tend = 0.0_wp
1824             CALL advec_s_bc( qr, 'qr' )
1825
1826          ENDIF
1827
1828!
1829!--       qr-tendency terms with no communication
1830          IF ( scalar_advec /= 'bc-scheme' )  THEN
1831             tend = 0.0_wp
1832             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1833                IF ( ws_scheme_sca )  THEN
1834                   CALL advec_s_ws( qr, 'qr' )
1835                ELSE
1836                   CALL advec_s_pw( qr )
1837                ENDIF
1838             ELSE
1839                CALL advec_s_up( qr )
1840             ENDIF
1841          ENDIF
1842
1843          CALL diffusion_s( qr,                                                &
1844                            surf_def_h(0)%qrsws, surf_def_h(1)%qrsws,          &
1845                            surf_def_h(2)%qrsws,                               &
1846                            surf_lsm_h%qrsws,    surf_usm_h%qrsws,             &
1847                            surf_def_v(0)%qrsws, surf_def_v(1)%qrsws,          &
1848                            surf_def_v(2)%qrsws, surf_def_v(3)%qrsws,          &
1849                            surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws,          &
1850                            surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws,          &
1851                            surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws,          &
1852                            surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws )
1853
1854!
1855!--       Prognostic equation for rain water content
1856          DO  i = nxl, nxr
1857             DO  j = nys, nyn
1858                DO  k = nzb+1, nzt
1859                   qr_p(k,j,i) = qr(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
1860                                                      tsc(3) * tqr_m(k,j,i) )  &
1861                                                    - tsc(5) * rdf_sc(k) *     &
1862                                                               qr(k,j,i)       &
1863                                             )                                 &
1864                                    * MERGE( 1.0_wp, 0.0_wp,                   &
1865                                             BTEST( wall_flags_0(k,j,i), 0 )   &
1866                                          )
1867                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
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                DO  i = nxl, nxr
1877                   DO  j = nys, nyn
1878                      DO  k = nzb+1, nzt
1879                         tqr_m(k,j,i) = tend(k,j,i)
1880                      ENDDO
1881                   ENDDO
1882                ENDDO
1883             ELSEIF ( intermediate_timestep_count < &
1884                      intermediate_timestep_count_max )  THEN
1885                DO  i = nxl, nxr
1886                   DO  j = nys, nyn
1887                      DO  k = nzb+1, nzt
1888                         tqr_m(k,j,i) =   -9.5625_wp * tend(k,j,i)             &
1889                                         + 5.3125_wp * tqr_m(k,j,i)
1890                      ENDDO
1891                   ENDDO
1892                ENDDO
1893             ENDIF
1894          ENDIF
1895
1896          CALL cpu_log( log_point(52), 'qr-equation', 'stop' )
1897          CALL cpu_log( log_point(53), 'nr-equation', 'start' )
1898
1899!
1900!--       Calculate prognostic equation for rain drop concentration
1901          sbt = tsc(2)
1902          IF ( scalar_advec == 'bc-scheme' )  THEN
1903
1904             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1905!
1906!--             Bott-Chlond scheme always uses Euler time step. Thus:
1907                sbt = 1.0_wp
1908             ENDIF
1909             tend = 0.0_wp
1910             CALL advec_s_bc( nr, 'nr' )
1911
1912          ENDIF
1913
1914!
1915!--       nr-tendency terms with no communication
1916          IF ( scalar_advec /= 'bc-scheme' )  THEN
1917             tend = 0.0_wp
1918             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1919                IF ( ws_scheme_sca )  THEN
1920                   CALL advec_s_ws( nr, 'nr' )
1921                ELSE
1922                   CALL advec_s_pw( nr )
1923                ENDIF
1924             ELSE
1925                CALL advec_s_up( nr )
1926             ENDIF
1927          ENDIF
1928
1929          CALL diffusion_s( nr,                                                &
1930                            surf_def_h(0)%nrsws, surf_def_h(1)%nrsws,          &
1931                            surf_def_h(2)%nrsws,                               &
1932                            surf_lsm_h%nrsws,    surf_usm_h%nrsws,             & 
1933                            surf_def_v(0)%nrsws, surf_def_v(1)%nrsws,          &
1934                            surf_def_v(2)%nrsws, surf_def_v(3)%nrsws,          &
1935                            surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws,          &
1936                            surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws,          &
1937                            surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws,          &
1938                            surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws )
1939
1940!
1941!--       Prognostic equation for rain drop concentration
1942          DO  i = nxl, nxr
1943             DO  j = nys, nyn
1944                DO  k = nzb+1, nzt
1945                   nr_p(k,j,i) = nr(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
1946                                                      tsc(3) * tnr_m(k,j,i) )  &
1947                                                    - tsc(5) * rdf_sc(k) *     &
1948                                                               nr(k,j,i)       &
1949                                             )                                 &
1950                                   * MERGE( 1.0_wp, 0.0_wp,                    &
1951                                             BTEST( wall_flags_0(k,j,i), 0 )   &
1952                                          )
1953                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
1954                ENDDO
1955             ENDDO
1956          ENDDO
1957
1958!
1959!--       Calculate tendencies for the next Runge-Kutta step
1960          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1961             IF ( intermediate_timestep_count == 1 )  THEN
1962                DO  i = nxl, nxr
1963                   DO  j = nys, nyn
1964                      DO  k = nzb+1, nzt
1965                         tnr_m(k,j,i) = tend(k,j,i)
1966                      ENDDO
1967                   ENDDO
1968                ENDDO
1969             ELSEIF ( intermediate_timestep_count < &
1970                      intermediate_timestep_count_max )  THEN
1971                DO  i = nxl, nxr
1972                   DO  j = nys, nyn
1973                      DO  k = nzb+1, nzt
1974                         tnr_m(k,j,i) =  -9.5625_wp * tend(k,j,i)             &
1975                                         + 5.3125_wp * tnr_m(k,j,i)
1976                      ENDDO
1977                   ENDDO
1978                ENDDO
1979             ENDIF
1980          ENDIF
1981
1982          CALL cpu_log( log_point(53), 'nr-equation', 'stop' )
1983
1984       ENDIF
1985
1986    ENDIF
1987!
1988!-- If required, compute prognostic equation for scalar
1989    IF ( passive_scalar )  THEN
1990
1991       CALL cpu_log( log_point(66), 's-equation', 'start' )
1992
1993!
1994!--    Scalar/q-tendency terms with communication
1995       sbt = tsc(2)
1996       IF ( scalar_advec == 'bc-scheme' )  THEN
1997
1998          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1999!
2000!--          Bott-Chlond scheme always uses Euler time step. Thus:
2001             sbt = 1.0_wp
2002          ENDIF
2003          tend = 0.0_wp
2004          CALL advec_s_bc( s, 's' )
2005
2006       ENDIF
2007
2008!
2009!--    Scalar/q-tendency terms with no communication
2010       IF ( scalar_advec /= 'bc-scheme' )  THEN
2011          tend = 0.0_wp
2012          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2013             IF ( ws_scheme_sca )  THEN
2014                CALL advec_s_ws( s, 's' )
2015             ELSE
2016                CALL advec_s_pw( s )
2017             ENDIF
2018          ELSE
2019             CALL advec_s_up( s )
2020          ENDIF
2021       ENDIF
2022
2023       CALL diffusion_s( s,                                                    &
2024                         surf_def_h(0)%ssws, surf_def_h(1)%ssws,               &
2025                         surf_def_h(2)%ssws,                                   &
2026                         surf_lsm_h%ssws,    surf_usm_h%ssws,                  &
2027                         surf_def_v(0)%ssws, surf_def_v(1)%ssws,               &
2028                         surf_def_v(2)%ssws, surf_def_v(3)%ssws,               &
2029                         surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,               &
2030                         surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,               &
2031                         surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,               &
2032                         surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
2033
2034!
2035!--    Sink or source of humidity due to canopy elements
2036       IF ( plant_canopy ) CALL pcm_tendency( 7 )
2037
2038!
2039!--    Large scale advection. Not implemented for scalars so far.
2040!        IF ( large_scale_forcing )  THEN
2041!           CALL ls_advec( simulated_time, 'q' )
2042!        ENDIF
2043
2044!
2045!--    Nudging. Not implemented for scalars so far.
2046!        IF ( nudging )  CALL nudge( simulated_time, 'q' )
2047
2048!
2049!--    If required compute influence of large-scale subsidence/ascent.
2050!--    Not implemented for scalars so far.
2051       IF ( large_scale_subsidence  .AND.                                      &
2052            .NOT. use_subsidence_tendencies  .AND.                             &
2053            .NOT. large_scale_forcing )  THEN
2054         CALL subsidence( tend, s, s_init, 3 )
2055       ENDIF
2056
2057       CALL user_actions( 's-tendency' )
2058
2059!
2060!--    Prognostic equation for total water content
2061       DO  i = nxl, nxr
2062          DO  j = nys, nyn
2063             DO  k = nzb+1, nzt
2064                s_p(k,j,i) = s(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
2065                                                 tsc(3) * ts_m(k,j,i) )        &
2066                                               - tsc(5) * rdf_sc(k) *          &
2067                                                    ( s(k,j,i) - s_init(k) )   &
2068                                        )                                      &
2069                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2070                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2071                                          )
2072                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
2073             ENDDO
2074          ENDDO
2075       ENDDO
2076
2077!
2078!--    Calculate tendencies for the next Runge-Kutta step
2079       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2080          IF ( intermediate_timestep_count == 1 )  THEN
2081             DO  i = nxl, nxr
2082                DO  j = nys, nyn
2083                   DO  k = nzb+1, nzt
2084                      ts_m(k,j,i) = tend(k,j,i)
2085                   ENDDO
2086                ENDDO
2087             ENDDO
2088          ELSEIF ( intermediate_timestep_count < &
2089                   intermediate_timestep_count_max )  THEN
2090             DO  i = nxl, nxr
2091                DO  j = nys, nyn
2092                   DO  k = nzb+1, nzt
2093                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
2094                                     + 5.3125_wp * ts_m(k,j,i)
2095                   ENDDO
2096                ENDDO
2097             ENDDO
2098          ENDIF
2099       ENDIF
2100
2101       CALL cpu_log( log_point(66), 's-equation', 'stop' )
2102
2103    ENDIF
2104!
2105!-- If required, compute prognostic equation for turbulent kinetic
2106!-- energy (TKE)
2107    IF ( .NOT. constant_diffusion )  THEN
2108
2109       CALL cpu_log( log_point(16), 'tke-equation', 'start' )
2110
2111       sbt = tsc(2)
2112       IF ( .NOT. use_upstream_for_tke )  THEN
2113          IF ( scalar_advec == 'bc-scheme' )  THEN
2114
2115             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2116!
2117!--             Bott-Chlond scheme always uses Euler time step. Thus:
2118                sbt = 1.0_wp
2119             ENDIF
2120             tend = 0.0_wp
2121             CALL advec_s_bc( e, 'e' )
2122
2123          ENDIF
2124       ENDIF
2125
2126!
2127!--    TKE-tendency terms with no communication
2128       IF ( scalar_advec /= 'bc-scheme'  .OR.  use_upstream_for_tke )  THEN
2129          IF ( use_upstream_for_tke )  THEN
2130             tend = 0.0_wp
2131             CALL advec_s_up( e )
2132          ELSE
2133             tend = 0.0_wp
2134             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2135                IF ( ws_scheme_sca )  THEN
2136                   CALL advec_s_ws( e, 'e' )
2137                ELSE
2138                   CALL advec_s_pw( e )
2139                ENDIF
2140             ELSE
2141                CALL advec_s_up( e )
2142             ENDIF
2143          ENDIF
2144       ENDIF
2145
2146       IF ( .NOT. humidity )  THEN
2147          IF ( ocean )  THEN
2148             CALL diffusion_e( prho, prho_reference )
2149          ELSE
2150             CALL diffusion_e( pt, pt_reference )
2151          ENDIF
2152       ELSE
2153          CALL diffusion_e( vpt, pt_reference )
2154       ENDIF
2155
2156       CALL production_e
2157
2158!
2159!--    Additional sink term for flows through plant canopies
2160       IF ( plant_canopy )  CALL pcm_tendency( 6 )
2161       CALL user_actions( 'e-tendency' )
2162
2163!
2164!--    Prognostic equation for TKE.
2165!--    Eliminate negative TKE values, which can occur due to numerical
2166!--    reasons in the course of the integration. In such cases the old TKE
2167!--    value is reduced by 90%.
2168       DO  i = nxl, nxr
2169          DO  j = nys, nyn
2170             DO  k = nzb+1, nzt
2171                e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
2172                                                 tsc(3) * te_m(k,j,i) )        &
2173                                        )                                      &
2174                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2175                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2176                                          )
2177                IF ( e_p(k,j,i) < 0.0_wp )  e_p(k,j,i) = 0.1_wp * e(k,j,i)
2178             ENDDO
2179          ENDDO
2180       ENDDO
2181
2182!
2183!--    Calculate tendencies for the next Runge-Kutta step
2184       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2185          IF ( intermediate_timestep_count == 1 )  THEN
2186             DO  i = nxl, nxr
2187                DO  j = nys, nyn
2188                   DO  k = nzb+1, nzt
2189                      te_m(k,j,i) = tend(k,j,i)
2190                   ENDDO
2191                ENDDO
2192             ENDDO
2193          ELSEIF ( intermediate_timestep_count < &
2194                   intermediate_timestep_count_max )  THEN
2195             DO  i = nxl, nxr
2196                DO  j = nys, nyn
2197                   DO  k = nzb+1, nzt
2198                      te_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
2199                                     + 5.3125_wp * te_m(k,j,i)
2200                   ENDDO
2201                ENDDO
2202             ENDDO
2203          ENDIF
2204       ENDIF
2205
2206       CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
2207
2208    ENDIF
2209
2210 END SUBROUTINE prognostic_equations_vector
2211
2212
2213 END MODULE prognostic_equations_mod
Note: See TracBrowser for help on using the repository browser.