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

Last change on this file since 2274 was 2261, checked in by raasch, 7 years ago

changes in mrun: unified cycle numbers for output files are used, paths and filenames are allowed to contain arbitrary numbers of dots, archive feature completely removed from the script, nech related parts completely removed, OpenMP bugfix in prognostic_equations

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