source: palm/trunk/SOURCE/check_parameters.f90 @ 2345

Last change on this file since 2345 was 2345, checked in by Giersch, 7 years ago

! Remove error message PA0156 and the conserve_volume_flow_mode option

  • Property svn:keywords set to Id
  • Property svn:mergeinfo set to False
    /palm/branches/forwind/SOURCE/check_parameters.f901564-1913
File size: 164.3 KB
Line 
1!> @file check_parameters.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: check_parameters.f90 2345 2017-08-09 11:50:30Z Giersch $
27! Remove error message PA0156 and the conserve_volume_flow_mode option
28! inflow_profile
29!
30! 2339 2017-08-07 13:55:26Z raasch
31! corrected timestamp in header
32!
33! 2338 2017-08-07 12:15:38Z gronemeier
34! Modularize 1D model
35!
36! 2329 2017-08-03 14:24:56Z knoop
37! Bugfix: index corrected for rho_air and rho_air_zw output
38!
39! 2320 2017-07-21 12:47:43Z suehring
40! Modularize large-scale forcing and nudging
41!
42! 2312 2017-07-14 20:26:51Z hoffmann
43! PA0349 and PA0420 removed.
44!
45! 2300 2017-06-29 13:31:14Z raasch
46! host-specific settings and checks removed
47!
48! 2292 2017-06-20 09:51:42Z schwenkel
49! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
50! includes two more prognostic equations for cloud drop concentration (nc)
51! and cloud water content (qc).
52!
53! 2274 2017-06-09 13:27:48Z Giersch
54! Changed error messages
55!
56! 2271 2017-06-09 12:34:55Z sward
57! roughness-length check altered
58! Error messages fixed
59!
60! 2270 2017-06-09 12:18:47Z maronga
61! Revised numbering (removed 2 timeseries)
62!
63! 2259 2017-06-08 09:09:11Z gronemeier
64! Implemented synthetic turbulence generator
65!
66! 2251 2017-06-06 15:10:46Z Giersch
67!
68! 2250 2017-06-06 15:07:12Z Giersch
69! Doxygen comment added
70!
71! 2248 2017-06-06 13:52:54Z sward
72! Error message fixed
73!
74! 2209 2017-04-19 09:34:46Z kanani
75! Check for plant canopy model output
76!
77! 2200 2017-04-11 11:37:51Z suehring
78! monotonic_adjustment removed
79!
80! 2178 2017-03-17 11:07:39Z hellstea
81! Index limits for perturbations are now set also in case of nested
82! boundary conditions
83!
84! 2169 2017-03-06 18:16:35Z suehring
85! Bugfix, move setting for topography grid convention to init_grid
86!
87! 2118 2017-01-17 16:38:49Z raasch
88! OpenACC related parts of code removed
89!
90! 2088 2016-12-19 16:30:25Z suehring
91! Bugfix in initial salinity profile
92!
93! 2084 2016-12-09 15:59:42Z knoop
94! Removed anelastic + multigrid error message
95!
96! 2052 2016-11-08 15:14:59Z gronemeier
97! Bugfix: remove setting of default value for recycling_width
98!
99! 2050 2016-11-08 15:00:55Z gronemeier
100! Implement turbulent outflow condition
101!
102! 2044 2016-11-02 16:44:25Z knoop
103! Added error code for anelastic approximation
104!
105! 2042 2016-11-02 13:47:31Z suehring
106! Additional checks for wall_heatflux, wall_humidityflux and wall_scalarflux.
107! Bugfix, check for constant_scalarflux.
108!
109! 2040 2016-10-26 16:58:09Z gronemeier
110! Removed check for statistic_regions > 9.
111!
112! 2037 2016-10-26 11:15:40Z knoop
113! Anelastic approximation implemented
114!
115! 2031 2016-10-21 15:11:58Z knoop
116! renamed variable rho to rho_ocean
117!
118! 2026 2016-10-18 10:27:02Z suehring
119! Bugfix, enable output of s*2
120!
121! 2024 2016-10-12 16:42:37Z kanani
122! Added missing CASE, error message and unit for ssws*,
123! increased number of possible output quantities to 500.
124!
125! 2011 2016-09-19 17:29:57Z kanani
126! Flag urban_surface is now defined in module control_parameters,
127! changed prefix for urban surface model output to "usm_",
128! added flag lsf_exception (inipar-Namelist parameter) to allow explicit
129! enabling of large scale forcing together with buildings on flat terrain,
130! introduced control parameter varnamelength for LEN of var.
131!
132! 2007 2016-08-24 15:47:17Z kanani
133! Added checks for the urban surface model,
134! increased counter in DO WHILE loop over data_output (for urban surface output)
135!
136! 2000 2016-08-20 18:09:15Z knoop
137! Forced header and separation lines into 80 columns
138!
139! 1994 2016-08-15 09:52:21Z suehring
140! Add missing check for cloud_physics and cloud_droplets
141!
142! 1992 2016-08-12 15:14:59Z suehring
143! New checks for top_scalarflux
144! Bugfixes concerning data output of profiles in case of passive_scalar
145!
146! 1984 2016-08-01 14:48:05Z suehring
147! Bugfix: checking for bottom and top boundary condition for humidity and scalars
148!
149! 1972 2016-07-26 07:52:02Z maronga
150! Removed check of lai* (done in land surface module)
151!
152! 1970 2016-07-18 14:27:12Z suehring
153! Bugfix, check for ibc_q_b and constant_waterflux in case of land-surface scheme
154!
155! 1962 2016-07-13 09:16:44Z suehring
156! typo removed
157!
158! 1960 2016-07-12 16:34:24Z suehring
159! Separate checks and calculations for humidity and passive scalar. Therefore,
160! a subroutine to calculate vertical gradients is introduced, in order  to reduce
161! redundance.
162! Additional check - large-scale forcing is not implemented for passive scalars
163! so far.
164!
165! 1931 2016-06-10 12:06:59Z suehring
166! Rename multigrid into multigrid_noopt and multigrid_fast into multigrid
167!
168! 1929 2016-06-09 16:25:25Z suehring
169! Bugfix in check for use_upstream_for_tke
170!
171! 1918 2016-05-27 14:35:57Z raasch
172! setting of a fixed reference state ('single_value') for ocean runs removed
173!
174! 1914 2016-05-26 14:44:07Z witha
175! Added checks for the wind turbine model
176!
177! 1841 2016-04-07 19:14:06Z raasch
178! redundant location message removed
179!
180! 1833 2016-04-07 14:23:03Z raasch
181! check of spectra quantities moved to spectra_mod
182!
183! 1829 2016-04-07 12:16:29Z maronga
184! Bugfix: output of user defined data required reset of the string 'unit'
185!
186! 1826 2016-04-07 12:01:39Z maronga
187! Moved checks for radiation model to the respective module.
188! Moved checks for plant canopy model to the respective module.
189! Bugfix for check of too large roughness_length
190!
191! 1824 2016-04-07 09:55:29Z gronemeier
192! Check if roughness_length < dz/2. Moved location_message(finished) to the end.
193!
194! 1822 2016-04-07 07:49:42Z hoffmann
195! PALM collision kernel deleted. Collision algorithms are checked for correct
196! spelling.
197!
198! Tails removed.
199! !
200! Checks for use_sgs_for_particles adopted for the use of droplets with
201! use_sgs_for_particles adopted.
202!
203! Unused variables removed.
204!
205! 1817 2016-04-06 15:44:20Z maronga
206! Moved checks for land_surface model to the respective module
207!
208! 1806 2016-04-05 18:55:35Z gronemeier
209! Check for recycling_yshift
210!
211! 1804 2016-04-05 16:30:18Z maronga
212! Removed code for parameter file check (__check)
213!
214! 1795 2016-03-18 15:00:53Z raasch
215! Bugfix: setting of initial scalar profile in ocean
216!
217! 1788 2016-03-10 11:01:04Z maronga
218! Added check for use of most_method = 'lookup' in combination with water
219! surface presribed in the land surface model. Added output of z0q.
220! Syntax layout improved.
221!
222! 1786 2016-03-08 05:49:27Z raasch
223! cpp-direktives for spectra removed, check of spectra level removed
224!
225! 1783 2016-03-06 18:36:17Z raasch
226! netcdf variables and module name changed,
227! check of netcdf precision removed (is done in the netcdf module)
228!
229! 1764 2016-02-28 12:45:19Z raasch
230! output of nest id in run description header,
231! bugfix: check of validity of lateral boundary conditions moved to parin
232!
233! 1762 2016-02-25 12:31:13Z hellstea
234! Introduction of nested domain feature
235!
236! 1745 2016-02-05 13:06:51Z gronemeier
237! Bugfix: check data output intervals to be /= 0.0 in case of parallel NetCDF4
238!
239! 1701 2015-11-02 07:43:04Z maronga
240! Bugfix: definition of rad_net timeseries was missing
241!
242! 1691 2015-10-26 16:17:44Z maronga
243! Added output of Obukhov length (ol) and radiative heating rates for RRTMG.
244! Added checks for use of radiation / lsm with topography.
245!
246! 1682 2015-10-07 23:56:08Z knoop
247! Code annotations made doxygen readable
248!
249! 1606 2015-06-29 10:43:37Z maronga
250! Added check for use of RRTMG without netCDF.
251!
252! 1587 2015-05-04 14:19:01Z maronga
253! Added check for set of albedos when using RRTMG
254!
255! 1585 2015-04-30 07:05:52Z maronga
256! Added support for RRTMG
257!
258! 1575 2015-03-27 09:56:27Z raasch
259! multigrid_fast added as allowed pressure solver
260!
261! 1573 2015-03-23 17:53:37Z suehring
262! Bugfix: check for advection schemes in case of non-cyclic boundary conditions
263!
264! 1557 2015-03-05 16:43:04Z suehring
265! Added checks for monotonic limiter
266!
267! 1555 2015-03-04 17:44:27Z maronga
268! Added output of r_a and r_s. Renumbering of LSM PA-messages.
269!
270! 1553 2015-03-03 17:33:54Z maronga
271! Removed check for missing soil temperature profile as default values were added.
272!
273! 1551 2015-03-03 14:18:16Z maronga
274! Added various checks for land surface and radiation model. In the course of this
275! action, the length of the variable var has to be increased
276!
277! 1504 2014-12-04 09:23:49Z maronga
278! Bugfix: check for large_scale forcing got lost.
279!
280! 1500 2014-12-03 17:42:41Z maronga
281! Boundary conditions changed to dirichlet for land surface model
282!
283! 1496 2014-12-02 17:25:50Z maronga
284! Added checks for the land surface model
285!
286! 1484 2014-10-21 10:53:05Z kanani
287! Changes due to new module structure of the plant canopy model:
288!   module plant_canopy_model_mod added,
289!   checks regarding usage of new method for leaf area density profile
290!   construction added,
291!   lad-profile construction moved to new subroutine init_plant_canopy within
292!   the module plant_canopy_model_mod,
293!   drag_coefficient renamed to canopy_drag_coeff.
294! Missing KIND-attribute for REAL constant added
295!
296! 1455 2014-08-29 10:47:47Z heinze
297! empty time records in volume, cross-section and masked data output prevented
298! in case of non-parallel netcdf-output in restart runs
299!
300! 1429 2014-07-15 12:53:45Z knoop
301! run_description_header exended to provide ensemble_member_nr if specified
302!
303! 1425 2014-07-05 10:57:53Z knoop
304! bugfix: perturbation domain modified for parallel random number generator
305!
306! 1402 2014-05-09 14:25:13Z raasch
307! location messages modified
308!
309! 1400 2014-05-09 14:03:54Z knoop
310! Check random generator extended by option random-parallel
311!
312! 1384 2014-05-02 14:31:06Z raasch
313! location messages added
314!
315! 1365 2014-04-22 15:03:56Z boeske
316! Usage of large scale forcing for pt and q enabled:
317! output for profiles of large scale advection (td_lsa_lpt, td_lsa_q),
318! large scale subsidence (td_sub_lpt, td_sub_q)
319! and nudging tendencies (td_nud_lpt, td_nud_q, td_nud_u and td_nud_v) added,
320! check if use_subsidence_tendencies is used correctly
321!
322! 1361 2014-04-16 15:17:48Z hoffmann
323! PA0363 removed
324! PA0362 changed
325!
326! 1359 2014-04-11 17:15:14Z hoffmann
327! Do not allow the execution of PALM with use_particle_tails, since particle
328! tails are currently not supported by our new particle structure.
329!
330! PA0084 not necessary for new particle structure
331!
332! 1353 2014-04-08 15:21:23Z heinze
333! REAL constants provided with KIND-attribute
334!
335! 1330 2014-03-24 17:29:32Z suehring
336! In case of SGS-particle velocity advection of TKE is also allowed with
337! dissipative 5th-order scheme.
338!
339! 1327 2014-03-21 11:00:16Z raasch
340! "baroclinicity" renamed "baroclinity", "ocean version" replaced by "ocean mode"
341! bugfix: duplicate error message 56 removed,
342! check of data_output_format and do3d_compress removed
343!
344! 1322 2014-03-20 16:38:49Z raasch
345! some REAL constants defined as wp-kind
346!
347! 1320 2014-03-20 08:40:49Z raasch
348! Kind-parameters added to all INTEGER and REAL declaration statements,
349! kinds are defined in new module kinds,
350! revision history before 2012 removed,
351! comment fields (!:) to be used for variable explanations added to
352! all variable declaration statements
353!
354! 1308 2014-03-13 14:58:42Z fricke
355! +netcdf_data_format_save
356! Calculate fixed number of output time levels for parallel netcdf output.
357! For masked data, parallel netcdf output is not tested so far, hence
358! netcdf_data_format is switched back to non-paralell output.
359!
360! 1299 2014-03-06 13:15:21Z heinze
361! enable usage of large_scale subsidence in combination with large_scale_forcing
362! output for profile of large scale vertical velocity w_subs added
363!
364! 1276 2014-01-15 13:40:41Z heinze
365! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
366!
367! 1241 2013-10-30 11:36:58Z heinze
368! output for profiles of ug and vg added
369! set surface_heatflux and surface_waterflux also in dependence on
370! large_scale_forcing
371! checks for nudging and large scale forcing from external file
372!
373! 1236 2013-09-27 07:21:13Z raasch
374! check number of spectra levels
375!
376! 1216 2013-08-26 09:31:42Z raasch
377! check for transpose_compute_overlap (temporary)
378!
379! 1214 2013-08-21 12:29:17Z kanani
380! additional check for simultaneous use of vertical grid stretching
381! and particle advection
382!
383! 1212 2013-08-15 08:46:27Z raasch
384! checks for poisfft_hybrid removed
385!
386! 1210 2013-08-14 10:58:20Z raasch
387! check for fftw
388!
389! 1179 2013-06-14 05:57:58Z raasch
390! checks and settings of buoyancy parameters and switches revised,
391! initial profile for rho_ocean added to hom (id=77)
392!
393! 1174 2013-05-31 10:28:08Z gryschka
394! Bugfix in computing initial profiles for ug, vg, lad, q in case of Atmosphere
395!
396! 1159 2013-05-21 11:58:22Z fricke
397! bc_lr/ns_dirneu/neudir removed
398!
399! 1115 2013-03-26 18:16:16Z hoffmann
400! unused variables removed
401! drizzle can be used without precipitation
402!
403! 1111 2013-03-08 23:54:10Z raasch
404! ibc_p_b = 2 removed
405!
406! 1103 2013-02-20 02:15:53Z raasch
407! Bugfix: turbulent inflow must not require cyclic fill in restart runs
408!
409! 1092 2013-02-02 11:24:22Z raasch
410! unused variables removed
411!
412! 1069 2012-11-28 16:18:43Z maronga
413! allow usage of topography in combination with cloud physics
414!
415! 1065 2012-11-22 17:42:36Z hoffmann
416! Bugfix: It is not allowed to use cloud_scheme = seifert_beheng without
417!         precipitation in order to save computational resources.
418!
419! 1060 2012-11-21 07:19:51Z raasch
420! additional check for parameter turbulent_inflow
421!
422! 1053 2012-11-13 17:11:03Z hoffmann
423! necessary changes for the new two-moment cloud physics scheme added:
424! - check cloud physics scheme (Kessler or Seifert and Beheng)
425! - plant_canopy is not allowed
426! - currently, only cache loop_optimization is allowed
427! - initial profiles of nr, qr
428! - boundary condition of nr, qr
429! - check output quantities (qr, nr, prr)
430!
431! 1036 2012-10-22 13:43:42Z raasch
432! code put under GPL (PALM 3.9)
433!
434! 1031/1034 2012-10-22 11:32:49Z raasch
435! check of netcdf4 parallel file support
436!
437! 1019 2012-09-28 06:46:45Z raasch
438! non-optimized version of prognostic_equations not allowed any more
439!
440! 1015 2012-09-27 09:23:24Z raasch
441! acc allowed for loop optimization,
442! checks for adjustment of mixing length to the Prandtl mixing length removed
443!
444! 1003 2012-09-14 14:35:53Z raasch
445! checks for cases with unequal subdomain sizes removed
446!
447! 1001 2012-09-13 14:08:46Z raasch
448! all actions concerning leapfrog- and upstream-spline-scheme removed
449!
450! 996 2012-09-07 10:41:47Z raasch
451! little reformatting
452
453! 978 2012-08-09 08:28:32Z fricke
454! setting of bc_lr/ns_dirneu/neudir
455! outflow damping layer removed
456! check for z0h*
457! check for pt_damping_width
458!
459! 964 2012-07-26 09:14:24Z raasch
460! check of old profil-parameters removed
461!
462! 940 2012-07-09 14:31:00Z raasch
463! checks for parameter neutral
464!
465! 924 2012-06-06 07:44:41Z maronga
466! Bugfix: preprocessor directives caused error during compilation
467!
468! 892 2012-05-02 13:51:44Z maronga
469! Bugfix for parameter file check ( excluding __netcdf4 )
470!
471! 866 2012-03-28 06:44:41Z raasch
472! use only 60% of the geostrophic wind as translation speed in case of Galilean
473! transformation and use_ug_for_galilei_tr = .T. in order to mimimize the
474! timestep
475!
476! 861 2012-03-26 14:18:34Z suehring
477! Check for topography and ws-scheme removed.
478! Check for loop_optimization = 'vector' and ws-scheme removed.
479!
480! 845 2012-03-07 10:23:05Z maronga
481! Bugfix: exclude __netcdf4 directive part from namelist file check compilation
482!
483! 828 2012-02-21 12:00:36Z raasch
484! check of collision_kernel extended
485!
486! 825 2012-02-19 03:03:44Z raasch
487! check for collision_kernel and curvature_solution_effects
488!
489! 809 2012-01-30 13:32:58Z maronga
490! Bugfix: replaced .AND. and .NOT. with && and ! in the preprocessor directives
491!
492! 807 2012-01-25 11:53:51Z maronga
493! New cpp directive "__check" implemented which is used by check_namelist_files
494!
495! Revision 1.1  1997/08/26 06:29:23  raasch
496! Initial revision
497!
498!
499! Description:
500! ------------
501!> Check control parameters and deduce further quantities.
502!------------------------------------------------------------------------------!
503 SUBROUTINE check_parameters
504
505
506    USE arrays_3d
507    USE cloud_parameters
508    USE constants
509    USE control_parameters
510    USE dvrp_variables
511    USE grid_variables
512    USE indices
513    USE land_surface_model_mod,                                                &
514        ONLY: lsm_check_data_output, lsm_check_data_output_pr,                 &
515              lsm_check_parameters
516
517    USE lsf_nudging_mod,                                                       &
518        ONLY:  lsf_nudging_check_parameters, lsf_nudging_check_data_output_pr, &
519               qsws_surf, shf_surf
520
521    USE kinds
522
523    USE model_1d_mod,                                                          &
524        ONLY:  damp_level_1d, damp_level_ind_1d
525
526    USE netcdf_interface,                                                      &
527        ONLY:  dopr_unit, do2d_unit, do3d_unit, netcdf_data_format,            &
528               netcdf_data_format_string, dots_unit, heatflux_output_unit,     &
529               waterflux_output_unit, momentumflux_output_unit
530
531    USE particle_attributes
532    USE pegrid
533    USE plant_canopy_model_mod,                                                &
534        ONLY:  pcm_check_data_output, pcm_check_parameters, plant_canopy
535
536    USE pmc_interface,                                                         &
537        ONLY:  cpl_id, nested_run
538
539    USE profil_parameter
540    USE radiation_model_mod,                                                   &
541        ONLY:  radiation, radiation_check_data_output,                         &
542               radiation_check_data_output_pr, radiation_check_parameters
543    USE spectra_mod,                                                           &
544        ONLY:  calculate_spectra, spectra_check_parameters
545    USE statistics
546    USE subsidence_mod
547    USE statistics
548    USE synthetic_turbulence_generator_mod,                                    &
549        ONLY:  stg_check_parameters
550    USE transpose_indices
551    USE urban_surface_mod,                                                     &
552        ONLY:  usm_check_data_output, usm_check_parameters
553    USE wind_turbine_model_mod,                                                &
554        ONLY:  wtm_check_parameters, wind_turbine
555
556
557    IMPLICIT NONE
558
559    CHARACTER (LEN=1)   ::  sq                       !<
560    CHARACTER (LEN=varnamelength)  ::  var           !<
561    CHARACTER (LEN=7)   ::  unit                     !<
562    CHARACTER (LEN=8)   ::  date                     !<
563    CHARACTER (LEN=10)  ::  time                     !<
564    CHARACTER (LEN=20)  ::  ensemble_string          !<
565    CHARACTER (LEN=15)  ::  nest_string              !<
566    CHARACTER (LEN=40)  ::  coupling_string          !<
567    CHARACTER (LEN=100) ::  action                   !<
568
569    INTEGER(iwp) ::  i                               !<
570    INTEGER(iwp) ::  ilen                            !<
571    INTEGER(iwp) ::  j                               !<
572    INTEGER(iwp) ::  k                               !<
573    INTEGER(iwp) ::  kk                              !<
574    INTEGER(iwp) ::  netcdf_data_format_save         !<
575    INTEGER(iwp) ::  position                        !<
576
577    LOGICAL     ::  found                            !<
578
579    REAL(wp)    ::  dum                              !<
580    REAL(wp)    ::  gradient                         !<
581    REAL(wp)    ::  remote = 0.0_wp                  !<
582    REAL(wp)    ::  simulation_time_since_reference  !<
583
584
585    CALL location_message( 'checking parameters', .FALSE. )
586
587!
588!-- Check for overlap combinations, which are not realized yet
589    IF ( transpose_compute_overlap )  THEN
590       IF ( numprocs == 1 )  STOP '+++ transpose-compute-overlap not implemented for single PE runs'
591    ENDIF
592
593!
594!-- Check the coupling mode
595!> @todo Check if any queries for other coupling modes (e.g. precursor_ocean) are missing
596    IF ( coupling_mode /= 'uncoupled'            .AND.  &
597         coupling_mode /= 'atmosphere_to_ocean'  .AND.  &
598         coupling_mode /= 'ocean_to_atmosphere' )  THEN
599       message_string = 'illegal coupling mode: ' // TRIM( coupling_mode )
600       CALL message( 'check_parameters', 'PA0002', 1, 2, 0, 6, 0 )
601    ENDIF
602
603!
604!-- Check dt_coupling, restart_time, dt_restart, end_time, dx, dy, nx and ny
605    IF ( coupling_mode /= 'uncoupled')  THEN
606
607       IF ( dt_coupling == 9999999.9_wp )  THEN
608          message_string = 'dt_coupling is not set but required for coup' //   &
609                           'ling mode "' //  TRIM( coupling_mode ) // '"'
610          CALL message( 'check_parameters', 'PA0003', 1, 2, 0, 6, 0 )
611       ENDIF
612
613#if defined( __parallel )
614
615!
616!--    NOTE: coupled runs have not been implemented in the check_namelist_files
617!--    program.
618!--    check_namelist_files will need the following information of the other
619!--    model (atmosphere/ocean).
620!       dt_coupling = remote
621!       dt_max = remote
622!       restart_time = remote
623!       dt_restart= remote
624!       simulation_time_since_reference = remote
625!       dx = remote
626
627
628       IF ( myid == 0 ) THEN
629          CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter,  &
630                         ierr )
631          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter,       &
632                         status, ierr )
633       ENDIF
634       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
635
636       IF ( dt_coupling /= remote )  THEN
637          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
638                 '": dt_coupling = ', dt_coupling, '& is not equal to ',       &
639                 'dt_coupling_remote = ', remote
640          CALL message( 'check_parameters', 'PA0004', 1, 2, 0, 6, 0 )
641       ENDIF
642       IF ( dt_coupling <= 0.0_wp )  THEN
643
644          IF ( myid == 0  ) THEN
645             CALL MPI_SEND( dt_max, 1, MPI_REAL, target_id, 19, comm_inter, ierr )
646             CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter,    &
647                            status, ierr )
648          ENDIF
649          CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
650
651          dt_coupling = MAX( dt_max, remote )
652          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
653                 '": dt_coupling <= 0.0 & is not allowed and is reset to ',    &
654                 'MAX(dt_max(A,O)) = ', dt_coupling
655          CALL message( 'check_parameters', 'PA0005', 0, 1, 0, 6, 0 )
656       ENDIF
657
658       IF ( myid == 0 ) THEN
659          CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, &
660                         ierr )
661          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter,       &
662                         status, ierr )
663       ENDIF
664       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
665
666       IF ( restart_time /= remote )  THEN
667          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
668                 '": restart_time = ', restart_time, '& is not equal to ',     &
669                 'restart_time_remote = ', remote
670          CALL message( 'check_parameters', 'PA0006', 1, 2, 0, 6, 0 )
671       ENDIF
672
673       IF ( myid == 0 ) THEN
674          CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter,   &
675                         ierr )
676          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter,       &
677                         status, ierr )
678       ENDIF
679       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
680
681       IF ( dt_restart /= remote )  THEN
682          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
683                 '": dt_restart = ', dt_restart, '& is not equal to ',         &
684                 'dt_restart_remote = ', remote
685          CALL message( 'check_parameters', 'PA0007', 1, 2, 0, 6, 0 )
686       ENDIF
687
688       simulation_time_since_reference = end_time - coupling_start_time
689
690       IF  ( myid == 0 ) THEN
691          CALL MPI_SEND( simulation_time_since_reference, 1, MPI_REAL, target_id, &
692                         14, comm_inter, ierr )
693          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter,       &
694                         status, ierr )
695       ENDIF
696       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
697
698       IF ( simulation_time_since_reference /= remote )  THEN
699          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
700                 '": simulation_time_since_reference = ',                      &
701                 simulation_time_since_reference, '& is not equal to ',        &
702                 'simulation_time_since_reference_remote = ', remote
703          CALL message( 'check_parameters', 'PA0008', 1, 2, 0, 6, 0 )
704       ENDIF
705
706       IF ( myid == 0 ) THEN
707          CALL MPI_SEND( dx, 1, MPI_REAL, target_id, 15, comm_inter, ierr )
708          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter,       &
709                                                             status, ierr )
710       ENDIF
711       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
712
713
714       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
715
716          IF ( dx < remote ) THEN
717             WRITE( message_string, * ) 'coupling mode "',                     &
718                   TRIM( coupling_mode ),                                      &
719           '": dx in Atmosphere is not equal to or not larger than dx in ocean'
720             CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 )
721          ENDIF
722
723          IF ( (nx_a+1)*dx /= (nx_o+1)*remote )  THEN
724             WRITE( message_string, * ) 'coupling mode "',                     &
725                    TRIM( coupling_mode ),                                     &
726             '": Domain size in x-direction is not equal in ocean and atmosphere'
727             CALL message( 'check_parameters', 'PA0010', 1, 2, 0, 6, 0 )
728          ENDIF
729
730       ENDIF
731
732       IF ( myid == 0) THEN
733          CALL MPI_SEND( dy, 1, MPI_REAL, target_id, 16, comm_inter, ierr )
734          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter,       &
735                         status, ierr )
736       ENDIF
737       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
738
739       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
740
741          IF ( dy < remote )  THEN
742             WRITE( message_string, * ) 'coupling mode "',                     &
743                    TRIM( coupling_mode ),                                     &
744                 '": dy in Atmosphere is not equal to or not larger than dy in ocean'
745             CALL message( 'check_parameters', 'PA0011', 1, 2, 0, 6, 0 )
746          ENDIF
747
748          IF ( (ny_a+1)*dy /= (ny_o+1)*remote )  THEN
749             WRITE( message_string, * ) 'coupling mode "',                     &
750                   TRIM( coupling_mode ),                                      &
751             '": Domain size in y-direction is not equal in ocean and atmosphere'
752             CALL message( 'check_parameters', 'PA0012', 1, 2, 0, 6, 0 )
753          ENDIF
754
755          IF ( MOD(nx_o+1,nx_a+1) /= 0 )  THEN
756             WRITE( message_string, * ) 'coupling mode "',                     &
757                   TRIM( coupling_mode ),                                      &
758             '": nx+1 in ocean is not divisible by nx+1 in', &
759             ' atmosphere without remainder'
760             CALL message( 'check_parameters', 'PA0339', 1, 2, 0, 6, 0 )
761          ENDIF
762
763          IF ( MOD(ny_o+1,ny_a+1) /= 0 )  THEN
764             WRITE( message_string, * ) 'coupling mode "',                     &
765                   TRIM( coupling_mode ),                                      &
766             '": ny+1 in ocean is not divisible by ny+1 in', &
767             ' atmosphere without remainder'
768
769             CALL message( 'check_parameters', 'PA0340', 1, 2, 0, 6, 0 )
770          ENDIF
771
772       ENDIF
773#else
774       WRITE( message_string, * ) 'coupling requires PALM to be called with',  &
775            ' ''mrun -K parallel'''
776       CALL message( 'check_parameters', 'PA0141', 1, 2, 0, 6, 0 )
777#endif
778    ENDIF
779
780#if defined( __parallel )
781!
782!-- Exchange via intercommunicator
783    IF ( coupling_mode == 'atmosphere_to_ocean' .AND. myid == 0 )  THEN
784       CALL MPI_SEND( humidity, 1, MPI_LOGICAL, target_id, 19, comm_inter,     &
785                      ierr )
786    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' .AND. myid == 0)  THEN
787       CALL MPI_RECV( humidity_remote, 1, MPI_LOGICAL, target_id, 19,          &
788                      comm_inter, status, ierr )
789    ENDIF
790    CALL MPI_BCAST( humidity_remote, 1, MPI_LOGICAL, 0, comm2d, ierr)
791
792#endif
793
794
795!
796!-- Generate the file header which is used as a header for most of PALM's
797!-- output files
798    CALL DATE_AND_TIME( date, time )
799    run_date = date(7:8)//'-'//date(5:6)//'-'//date(3:4)
800    run_time = time(1:2)//':'//time(3:4)//':'//time(5:6)
801    IF ( coupling_mode == 'uncoupled' )  THEN
802       coupling_string = ''
803    ELSEIF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
804       coupling_string = ' coupled (atmosphere)'
805    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
806       coupling_string = ' coupled (ocean)'
807    ENDIF
808    IF ( ensemble_member_nr /= 0 )  THEN
809       WRITE( ensemble_string, '(2X,A,I2.2)' )  'en-no: ', ensemble_member_nr
810    ELSE
811       ensemble_string = ''
812    ENDIF
813    IF ( nested_run )  THEN
814       WRITE( nest_string, '(2X,A,I2.2)' )  'nest-id: ', cpl_id
815    ELSE
816       nest_string = ''
817    ENDIF
818
819    WRITE ( run_description_header,                                            &
820            '(A,2X,A,2X,A,A,A,I2.2,A,A,A,2X,A,A,2X,A,1X,A)' )                  &
821          TRIM( version ), TRIM( revision ), 'run: ',                          &
822          TRIM( run_identifier ), '.', runnr, TRIM( coupling_string ),         &
823          TRIM( nest_string ), TRIM( ensemble_string), 'host: ', TRIM( host ), &
824          run_date, run_time
825
826!
827!-- Check the general loop optimization method
828    SELECT CASE ( TRIM( loop_optimization ) )
829
830       CASE ( 'cache', 'vector' )
831          CONTINUE
832
833       CASE DEFAULT
834          message_string = 'illegal value given for loop_optimization: "' //   &
835                           TRIM( loop_optimization ) // '"'
836          CALL message( 'check_parameters', 'PA0013', 1, 2, 0, 6, 0 )
837
838    END SELECT
839
840!
841!-- Check if vertical grid stretching is used together with particles
842    IF ( dz_stretch_level < 100000.0_wp .AND. particle_advection )  THEN
843       message_string = 'Vertical grid stretching is not allowed together ' // &
844                        'with particle advection.'
845       CALL message( 'check_parameters', 'PA0017', 1, 2, 0, 6, 0 )
846    ENDIF
847
848!
849!-- Check topography setting (check for illegal parameter combinations)
850    IF ( topography /= 'flat' )  THEN
851       action = ' '
852       IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme'      &
853          )  THEN
854          WRITE( action, '(A,A)' )  'scalar_advec = ', scalar_advec
855       ENDIF
856       IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' )&
857       THEN
858          WRITE( action, '(A,A)' )  'momentum_advec = ', momentum_advec
859       ENDIF
860       IF ( psolver == 'sor' )  THEN
861          WRITE( action, '(A,A)' )  'psolver = ', psolver
862       ENDIF
863       IF ( sloping_surface )  THEN
864          WRITE( action, '(A)' )  'sloping surface = .TRUE.'
865       ENDIF
866       IF ( galilei_transformation )  THEN
867          WRITE( action, '(A)' )  'galilei_transformation = .TRUE.'
868       ENDIF
869       IF ( cloud_physics )  THEN
870          WRITE( action, '(A)' )  'cloud_physics = .TRUE.'
871       ENDIF
872       IF ( cloud_droplets )  THEN
873          WRITE( action, '(A)' )  'cloud_droplets = .TRUE.'
874       ENDIF
875       IF ( .NOT. constant_flux_layer )  THEN
876          WRITE( action, '(A)' )  'constant_flux_layer = .FALSE.'
877       ENDIF
878       IF ( action /= ' ' )  THEN
879          message_string = 'a non-flat topography does not allow ' //          &
880                           TRIM( action )
881          CALL message( 'check_parameters', 'PA0014', 1, 2, 0, 6, 0 )
882       ENDIF
883
884    ENDIF
885
886!
887!-- Check approximation
888    IF ( TRIM( approximation ) /= 'boussinesq'   .AND.   &
889         TRIM( approximation ) /= 'anelastic' )  THEN
890       message_string = 'unknown approximation: approximation = "' //    &
891                        TRIM( approximation ) // '"'
892       CALL message( 'check_parameters', 'PA0446', 1, 2, 0, 6, 0 )
893    ENDIF
894
895!
896!-- Check approximation requirements
897    IF ( TRIM( approximation ) == 'anelastic'   .AND.   &
898         TRIM( momentum_advec ) /= 'ws-scheme' )  THEN
899       message_string = 'Anelastic approximation requires: ' //                &
900                        'momentum_advec = "ws-scheme"'
901       CALL message( 'check_parameters', 'PA0447', 1, 2, 0, 6, 0 )
902    ENDIF
903    IF ( TRIM( approximation ) == 'anelastic'   .AND.   &
904         TRIM( psolver ) == 'multigrid' )  THEN
905       message_string = 'Anelastic approximation currently only supports: ' // &
906                        'psolver = "poisfft", ' // &
907                        'psolver = "sor" and ' // &
908                        'psolver = "multigrid_noopt"'
909       CALL message( 'check_parameters', 'PA0448', 1, 2, 0, 6, 0 )
910    ENDIF
911    IF ( TRIM( approximation ) == 'anelastic'   .AND.   &
912         conserve_volume_flow )  THEN
913       message_string = 'Anelastic approximation is not allowed with:' //      &
914                        'conserve_volume_flow = .TRUE.'
915       CALL message( 'check_parameters', 'PA0449', 1, 2, 0, 6, 0 )
916    ENDIF
917
918!
919!-- Check flux input mode
920    IF ( TRIM( flux_input_mode ) /= 'dynamic'    .AND.   &
921         TRIM( flux_input_mode ) /= 'kinematic'  .AND.   &
922         TRIM( flux_input_mode ) /= 'approximation-specific' )  THEN
923       message_string = 'unknown flux input mode: flux_input_mode = "' //      &
924                        TRIM( flux_input_mode ) // '"'
925       CALL message( 'check_parameters', 'PA0450', 1, 2, 0, 6, 0 )
926    ENDIF
927!-- Set flux input mode according to approximation if applicable
928    IF ( TRIM( flux_input_mode ) == 'approximation-specific' )  THEN
929       IF ( TRIM( approximation ) == 'anelastic' )  THEN
930          flux_input_mode = 'dynamic'
931       ELSEIF ( TRIM( approximation ) == 'boussinesq' )  THEN
932          flux_input_mode = 'kinematic'
933       ENDIF
934    ENDIF
935
936!
937!-- Check flux output mode
938    IF ( TRIM( flux_output_mode ) /= 'dynamic'    .AND.   &
939         TRIM( flux_output_mode ) /= 'kinematic'  .AND.   &
940         TRIM( flux_output_mode ) /= 'approximation-specific' )  THEN
941       message_string = 'unknown flux output mode: flux_output_mode = "' //    &
942                        TRIM( flux_output_mode ) // '"'
943       CALL message( 'check_parameters', 'PA0451', 1, 2, 0, 6, 0 )
944    ENDIF
945!-- Set flux output mode according to approximation if applicable
946    IF ( TRIM( flux_output_mode ) == 'approximation-specific' )  THEN
947       IF ( TRIM( approximation ) == 'anelastic' )  THEN
948          flux_output_mode = 'dynamic'
949       ELSEIF ( TRIM( approximation ) == 'boussinesq' )  THEN
950          flux_output_mode = 'kinematic'
951       ENDIF
952    ENDIF
953
954
955!
956!-- When the land surface model is used, the flux output must be dynamic.
957    IF ( land_surface )  THEN
958       flux_output_mode = 'dynamic'
959    ENDIF
960
961!
962!-- set the flux output units according to flux_output_mode
963    IF ( TRIM( flux_output_mode ) == 'kinematic' ) THEN
964        heatflux_output_unit              = 'K m/s'
965        waterflux_output_unit             = 'kg/kg m/s'
966        momentumflux_output_unit          = 'm2/s2'
967    ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' ) THEN
968        heatflux_output_unit              = 'W/m2'
969        waterflux_output_unit             = 'W/m2'
970        momentumflux_output_unit          = 'N/m2'
971    ENDIF
972
973!-- set time series output units for fluxes
974    dots_unit(14:16) = heatflux_output_unit
975    dots_unit(21)    = waterflux_output_unit
976    dots_unit(19:20) = momentumflux_output_unit
977
978!
979!-- Check ocean setting
980    IF ( TRIM( coupling_mode ) == 'uncoupled'  .AND.                           &
981         TRIM( coupling_char ) == '_O' .AND.                                   &
982         .NOT. ocean )  THEN
983
984!
985!--    Check whether an (uncoupled) atmospheric run has been declared as an
986!--    ocean run (this setting is done via mrun-option -y)
987       message_string = 'ocean = .F. does not allow coupling_char = "' //      &
988                        TRIM( coupling_char ) // '" set by mrun-option "-y"'
989       CALL message( 'check_parameters', 'PA0317', 1, 2, 0, 6, 0 )
990
991    ENDIF
992!
993!-- Check cloud scheme
994    IF ( cloud_scheme == 'saturation_adjust' )  THEN
995       microphysics_sat_adjust = .TRUE.
996       microphysics_seifert    = .FALSE.
997       microphysics_kessler    = .FALSE.
998       precipitation           = .FALSE.
999    ELSEIF ( cloud_scheme == 'seifert_beheng' )  THEN
1000       microphysics_sat_adjust = .FALSE.
1001       microphysics_seifert    = .TRUE.
1002       microphysics_kessler    = .FALSE.
1003       microphysics_morrison  = .FALSE.
1004       precipitation           = .TRUE.
1005    ELSEIF ( cloud_scheme == 'kessler' )  THEN
1006       microphysics_sat_adjust = .FALSE.
1007       microphysics_seifert    = .FALSE.
1008       microphysics_kessler    = .TRUE.
1009       microphysics_morrison   = .FALSE.
1010       precipitation           = .TRUE.
1011    ELSEIF ( cloud_scheme == 'morrison' )  THEN
1012       microphysics_sat_adjust = .FALSE.
1013       microphysics_seifert    = .TRUE.
1014       microphysics_kessler    = .FALSE.
1015       microphysics_morrison   = .TRUE.
1016       precipitation           = .TRUE.
1017    ELSE
1018       message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // &
1019                        TRIM( cloud_scheme ) // '"'
1020       CALL message( 'check_parameters', 'PA0357', 1, 2, 0, 6, 0 )
1021    ENDIF
1022!
1023!-- Check whether there are any illegal values
1024!-- Pressure solver:
1025    IF ( psolver /= 'poisfft'  .AND.  psolver /= 'sor'  .AND.                  &
1026         psolver /= 'multigrid'  .AND.  psolver /= 'multigrid_noopt' )  THEN
1027       message_string = 'unknown solver for perturbation pressure: psolver' // &
1028                        ' = "' // TRIM( psolver ) // '"'
1029       CALL message( 'check_parameters', 'PA0016', 1, 2, 0, 6, 0 )
1030    ENDIF
1031
1032    IF ( psolver(1:9) == 'multigrid' )  THEN
1033       IF ( cycle_mg == 'w' )  THEN
1034          gamma_mg = 2
1035       ELSEIF ( cycle_mg == 'v' )  THEN
1036          gamma_mg = 1
1037       ELSE
1038          message_string = 'unknown multigrid cycle: cycle_mg = "' //          &
1039                           TRIM( cycle_mg ) // '"'
1040          CALL message( 'check_parameters', 'PA0020', 1, 2, 0, 6, 0 )
1041       ENDIF
1042    ENDIF
1043
1044    IF ( fft_method /= 'singleton-algorithm'  .AND.                            &
1045         fft_method /= 'temperton-algorithm'  .AND.                            &
1046         fft_method /= 'fftw'                 .AND.                            &
1047         fft_method /= 'system-specific' )  THEN
1048       message_string = 'unknown fft-algorithm: fft_method = "' //             &
1049                        TRIM( fft_method ) // '"'
1050       CALL message( 'check_parameters', 'PA0021', 1, 2, 0, 6, 0 )
1051    ENDIF
1052
1053    IF( momentum_advec == 'ws-scheme' .AND.                                    &
1054        .NOT. call_psolver_at_all_substeps  ) THEN
1055        message_string = 'psolver must be called at each RK3 substep when "'// &
1056                      TRIM(momentum_advec) // ' "is used for momentum_advec'
1057        CALL message( 'check_parameters', 'PA0344', 1, 2, 0, 6, 0 )
1058    END IF
1059!
1060!-- Advection schemes:
1061    IF ( momentum_advec /= 'pw-scheme'  .AND.  momentum_advec /= 'ws-scheme' ) &
1062    THEN
1063       message_string = 'unknown advection scheme: momentum_advec = "' //      &
1064                        TRIM( momentum_advec ) // '"'
1065       CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 )
1066    ENDIF
1067    IF ( ( momentum_advec == 'ws-scheme' .OR.  scalar_advec == 'ws-scheme' )   &
1068           .AND. ( timestep_scheme == 'euler' .OR.                             &
1069                   timestep_scheme == 'runge-kutta-2' ) )                      &
1070    THEN
1071       message_string = 'momentum_advec or scalar_advec = "'                   &
1072         // TRIM( momentum_advec ) // '" is not allowed with timestep_scheme = "' // &
1073         TRIM( timestep_scheme ) // '"'
1074       CALL message( 'check_parameters', 'PA0023', 1, 2, 0, 6, 0 )
1075    ENDIF
1076    IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'ws-scheme' .AND. &
1077         scalar_advec /= 'bc-scheme' )                                         &
1078    THEN
1079       message_string = 'unknown advection scheme: scalar_advec = "' //        &
1080                        TRIM( scalar_advec ) // '"'
1081       CALL message( 'check_parameters', 'PA0024', 1, 2, 0, 6, 0 )
1082    ENDIF
1083    IF ( scalar_advec == 'bc-scheme'  .AND.  loop_optimization == 'cache' )    &
1084    THEN
1085       message_string = 'advection_scheme scalar_advec = "'                    &
1086         // TRIM( scalar_advec ) // '" not implemented for & loop_optimization = "' // &
1087         TRIM( loop_optimization ) // '"'
1088       CALL message( 'check_parameters', 'PA0026', 1, 2, 0, 6, 0 )
1089    ENDIF
1090
1091    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets  .AND.             &
1092         .NOT. use_upstream_for_tke  .AND.                                     &
1093         scalar_advec /= 'ws-scheme'                                           &
1094       )  THEN
1095       use_upstream_for_tke = .TRUE.
1096       message_string = 'use_upstream_for_tke is set to .TRUE. because ' //    &
1097                        'use_sgs_for_particles = .TRUE. '          //          &
1098                        'and scalar_advec /= ws-scheme'
1099       CALL message( 'check_parameters', 'PA0025', 0, 1, 0, 6, 0 )
1100    ENDIF
1101
1102!
1103!-- Set LOGICAL switches to enhance performance
1104    IF ( momentum_advec == 'ws-scheme' )  ws_scheme_mom = .TRUE.
1105    IF ( scalar_advec   == 'ws-scheme' )  ws_scheme_sca = .TRUE.
1106
1107
1108!
1109!-- Timestep schemes:
1110    SELECT CASE ( TRIM( timestep_scheme ) )
1111
1112       CASE ( 'euler' )
1113          intermediate_timestep_count_max = 1
1114
1115       CASE ( 'runge-kutta-2' )
1116          intermediate_timestep_count_max = 2
1117
1118       CASE ( 'runge-kutta-3' )
1119          intermediate_timestep_count_max = 3
1120
1121       CASE DEFAULT
1122          message_string = 'unknown timestep scheme: timestep_scheme = "' //   &
1123                           TRIM( timestep_scheme ) // '"'
1124          CALL message( 'check_parameters', 'PA0027', 1, 2, 0, 6, 0 )
1125
1126    END SELECT
1127
1128    IF ( (momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme')   &
1129         .AND. timestep_scheme(1:5) == 'runge' ) THEN
1130       message_string = 'momentum advection scheme "' // &
1131                        TRIM( momentum_advec ) // '" & does not work with ' // &
1132                        'timestep_scheme "' // TRIM( timestep_scheme ) // '"'
1133       CALL message( 'check_parameters', 'PA0029', 1, 2, 0, 6, 0 )
1134    ENDIF
1135!
1136!-- Check for proper settings for microphysics
1137    IF ( cloud_physics  .AND.  cloud_droplets )  THEN
1138       message_string = 'cloud_physics = .TRUE. is not allowed with ' //  &
1139                        'cloud_droplets = .TRUE.'
1140       CALL message( 'check_parameters', 'PA0442', 1, 2, 0, 6, 0 )
1141    ENDIF
1142
1143!
1144!-- Collision kernels:
1145    SELECT CASE ( TRIM( collision_kernel ) )
1146
1147       CASE ( 'hall', 'hall_fast' )
1148          hall_kernel = .TRUE.
1149
1150       CASE ( 'wang', 'wang_fast' )
1151          wang_kernel = .TRUE.
1152
1153       CASE ( 'none' )
1154
1155
1156       CASE DEFAULT
1157          message_string = 'unknown collision kernel: collision_kernel = "' // &
1158                           TRIM( collision_kernel ) // '"'
1159          CALL message( 'check_parameters', 'PA0350', 1, 2, 0, 6, 0 )
1160
1161    END SELECT
1162    IF ( collision_kernel(6:9) == 'fast' )  use_kernel_tables = .TRUE.
1163
1164    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
1165         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
1166!
1167!--    No restart run: several initialising actions are possible
1168       action = initializing_actions
1169       DO  WHILE ( TRIM( action ) /= '' )
1170          position = INDEX( action, ' ' )
1171          SELECT CASE ( action(1:position-1) )
1172
1173             CASE ( 'set_constant_profiles', 'set_1d-model_profiles',          &
1174                    'by_user', 'initialize_vortex', 'initialize_ptanom' )
1175                action = action(position+1:)
1176
1177             CASE DEFAULT
1178                message_string = 'initializing_action = "' //                  &
1179                                 TRIM( action ) // '" unkown or not allowed'
1180                CALL message( 'check_parameters', 'PA0030', 1, 2, 0, 6, 0 )
1181
1182          END SELECT
1183       ENDDO
1184    ENDIF
1185
1186    IF ( TRIM( initializing_actions ) == 'initialize_vortex'  .AND.            &
1187         conserve_volume_flow ) THEN
1188         message_string = 'initializing_actions = "initialize_vortex"' //      &
1189                        ' is not allowed with conserve_volume_flow = .T.'
1190       CALL message( 'check_parameters', 'PA0343', 1, 2, 0, 6, 0 )
1191    ENDIF
1192
1193
1194    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND.    &
1195         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
1196       message_string = 'initializing_actions = "set_constant_profiles"' //    &
1197                        ' and "set_1d-model_profiles" are not allowed ' //     &
1198                        'simultaneously'
1199       CALL message( 'check_parameters', 'PA0031', 1, 2, 0, 6, 0 )
1200    ENDIF
1201
1202    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND.    &
1203         INDEX( initializing_actions, 'by_user' ) /= 0 )  THEN
1204       message_string = 'initializing_actions = "set_constant_profiles"' //    &
1205                        ' and "by_user" are not allowed simultaneously'
1206       CALL message( 'check_parameters', 'PA0032', 1, 2, 0, 6, 0 )
1207    ENDIF
1208
1209    IF ( INDEX( initializing_actions, 'by_user' ) /= 0  .AND.                  &
1210         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
1211       message_string = 'initializing_actions = "by_user" and ' //             &
1212                        '"set_1d-model_profiles" are not allowed simultaneously'
1213       CALL message( 'check_parameters', 'PA0033', 1, 2, 0, 6, 0 )
1214    ENDIF
1215
1216    IF ( cloud_physics  .AND.  .NOT.  humidity )  THEN
1217       WRITE( message_string, * ) 'cloud_physics = ', cloud_physics, ' is ',   &
1218              'not allowed with humidity = ', humidity
1219       CALL message( 'check_parameters', 'PA0034', 1, 2, 0, 6, 0 )
1220    ENDIF
1221
1222    IF ( humidity  .AND.  sloping_surface )  THEN
1223       message_string = 'humidity = .TRUE. and sloping_surface = .TRUE. ' //   &
1224                        'are not allowed simultaneously'
1225       CALL message( 'check_parameters', 'PA0036', 1, 2, 0, 6, 0 )
1226    ENDIF
1227
1228!
1229!-- When synthetic turbulence generator is used, peform addtional checks
1230    IF ( synthetic_turbulence_generator )  CALL stg_check_parameters
1231!
1232!-- When plant canopy model is used, peform addtional checks
1233    IF ( plant_canopy )  CALL pcm_check_parameters
1234
1235!
1236!-- Additional checks for spectra
1237    IF ( calculate_spectra )  CALL spectra_check_parameters
1238!
1239!-- When land surface model is used, perform additional checks
1240    IF ( land_surface )  CALL lsm_check_parameters
1241
1242!
1243!-- When urban surface model is used, perform additional checks
1244    IF ( urban_surface )  CALL usm_check_parameters
1245
1246!
1247!-- If wind turbine model is used, perform additional checks
1248    IF ( wind_turbine )  CALL wtm_check_parameters
1249!
1250!-- When radiation model is used, peform addtional checks
1251    IF ( radiation )  CALL radiation_check_parameters
1252!
1253!-- When large-scale forcing or nudging is used, peform addtional checks
1254    IF ( large_scale_forcing  .OR.  nudging )  CALL lsf_nudging_check_parameters
1255
1256
1257
1258    IF ( .NOT. ( loop_optimization == 'cache'  .OR.                            &
1259                 loop_optimization == 'vector' )                               &
1260         .AND.  cloud_physics  .AND.  microphysics_seifert )  THEN
1261       message_string = 'cloud_scheme = seifert_beheng requires ' //           &
1262                        'loop_optimization = "cache" or "vector"'
1263       CALL message( 'check_parameters', 'PA0362', 1, 2, 0, 6, 0 )
1264    ENDIF
1265
1266!
1267!-- In case of no model continuation run, check initialising parameters and
1268!-- deduce further quantities
1269    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1270
1271!
1272!--    Initial profiles for 1D and 3D model, respectively (u,v further below)
1273       pt_init = pt_surface
1274       IF ( humidity       )  q_init  = q_surface
1275       IF ( ocean          )  sa_init = sa_surface
1276       IF ( passive_scalar )  s_init  = s_surface
1277
1278!
1279!--
1280!--    If required, compute initial profile of the geostrophic wind
1281!--    (component ug)
1282       i = 1
1283       gradient = 0.0_wp
1284
1285       IF (  .NOT.  ocean )  THEN
1286
1287          ug_vertical_gradient_level_ind(1) = 0
1288          ug(0) = ug_surface
1289          DO  k = 1, nzt+1
1290             IF ( i < 11 )  THEN
1291                IF ( ug_vertical_gradient_level(i) < zu(k)  .AND.              &
1292                     ug_vertical_gradient_level(i) >= 0.0_wp )  THEN
1293                   gradient = ug_vertical_gradient(i) / 100.0_wp
1294                   ug_vertical_gradient_level_ind(i) = k - 1
1295                   i = i + 1
1296                ENDIF
1297             ENDIF
1298             IF ( gradient /= 0.0_wp )  THEN
1299                IF ( k /= 1 )  THEN
1300                   ug(k) = ug(k-1) + dzu(k) * gradient
1301                ELSE
1302                   ug(k) = ug_surface + dzu(k) * gradient
1303                ENDIF
1304             ELSE
1305                ug(k) = ug(k-1)
1306             ENDIF
1307          ENDDO
1308
1309       ELSE
1310
1311          ug_vertical_gradient_level_ind(1) = nzt+1
1312          ug(nzt+1) = ug_surface
1313          DO  k = nzt, nzb, -1
1314             IF ( i < 11 )  THEN
1315                IF ( ug_vertical_gradient_level(i) > zu(k)  .AND.              &
1316                     ug_vertical_gradient_level(i) <= 0.0_wp )  THEN
1317                   gradient = ug_vertical_gradient(i) / 100.0_wp
1318                   ug_vertical_gradient_level_ind(i) = k + 1
1319                   i = i + 1
1320                ENDIF
1321             ENDIF
1322             IF ( gradient /= 0.0_wp )  THEN
1323                IF ( k /= nzt )  THEN
1324                   ug(k) = ug(k+1) - dzu(k+1) * gradient
1325                ELSE
1326                   ug(k)   = ug_surface - 0.5_wp * dzu(k+1) * gradient
1327                   ug(k+1) = ug_surface + 0.5_wp * dzu(k+1) * gradient
1328                ENDIF
1329             ELSE
1330                ug(k) = ug(k+1)
1331             ENDIF
1332          ENDDO
1333
1334       ENDIF
1335
1336!
1337!--    In case of no given gradients for ug, choose a zero gradient
1338       IF ( ug_vertical_gradient_level(1) == -9999999.9_wp )  THEN
1339          ug_vertical_gradient_level(1) = 0.0_wp
1340       ENDIF
1341
1342!
1343!--
1344!--    If required, compute initial profile of the geostrophic wind
1345!--    (component vg)
1346       i = 1
1347       gradient = 0.0_wp
1348
1349       IF (  .NOT.  ocean )  THEN
1350
1351          vg_vertical_gradient_level_ind(1) = 0
1352          vg(0) = vg_surface
1353          DO  k = 1, nzt+1
1354             IF ( i < 11 )  THEN
1355                IF ( vg_vertical_gradient_level(i) < zu(k)  .AND.              &
1356                     vg_vertical_gradient_level(i) >= 0.0_wp )  THEN
1357                   gradient = vg_vertical_gradient(i) / 100.0_wp
1358                   vg_vertical_gradient_level_ind(i) = k - 1
1359                   i = i + 1
1360                ENDIF
1361             ENDIF
1362             IF ( gradient /= 0.0_wp )  THEN
1363                IF ( k /= 1 )  THEN
1364                   vg(k) = vg(k-1) + dzu(k) * gradient
1365                ELSE
1366                   vg(k) = vg_surface + dzu(k) * gradient
1367                ENDIF
1368             ELSE
1369                vg(k) = vg(k-1)
1370             ENDIF
1371          ENDDO
1372
1373       ELSE
1374
1375          vg_vertical_gradient_level_ind(1) = nzt+1
1376          vg(nzt+1) = vg_surface
1377          DO  k = nzt, nzb, -1
1378             IF ( i < 11 )  THEN
1379                IF ( vg_vertical_gradient_level(i) > zu(k)  .AND.              &
1380                     vg_vertical_gradient_level(i) <= 0.0_wp )  THEN
1381                   gradient = vg_vertical_gradient(i) / 100.0_wp
1382                   vg_vertical_gradient_level_ind(i) = k + 1
1383                   i = i + 1
1384                ENDIF
1385             ENDIF
1386             IF ( gradient /= 0.0_wp )  THEN
1387                IF ( k /= nzt )  THEN
1388                   vg(k) = vg(k+1) - dzu(k+1) * gradient
1389                ELSE
1390                   vg(k)   = vg_surface - 0.5_wp * dzu(k+1) * gradient
1391                   vg(k+1) = vg_surface + 0.5_wp * dzu(k+1) * gradient
1392                ENDIF
1393             ELSE
1394                vg(k) = vg(k+1)
1395             ENDIF
1396          ENDDO
1397
1398       ENDIF
1399
1400!
1401!--    In case of no given gradients for vg, choose a zero gradient
1402       IF ( vg_vertical_gradient_level(1) == -9999999.9_wp )  THEN
1403          vg_vertical_gradient_level(1) = 0.0_wp
1404       ENDIF
1405
1406!
1407!--    Let the initial wind profiles be the calculated ug/vg profiles or
1408!--    interpolate them from wind profile data (if given)
1409       IF ( u_profile(1) == 9999999.9_wp  .AND.  v_profile(1) == 9999999.9_wp )  THEN
1410
1411          u_init = ug
1412          v_init = vg
1413
1414       ELSEIF ( u_profile(1) == 0.0_wp  .AND.  v_profile(1) == 0.0_wp )  THEN
1415
1416          IF ( uv_heights(1) /= 0.0_wp )  THEN
1417             message_string = 'uv_heights(1) must be 0.0'
1418             CALL message( 'check_parameters', 'PA0345', 1, 2, 0, 6, 0 )
1419          ENDIF
1420
1421          use_prescribed_profile_data = .TRUE.
1422
1423          kk = 1
1424          u_init(0) = 0.0_wp
1425          v_init(0) = 0.0_wp
1426
1427          DO  k = 1, nz+1
1428
1429             IF ( kk < 100 )  THEN
1430                DO  WHILE ( uv_heights(kk+1) <= zu(k) )
1431                   kk = kk + 1
1432                   IF ( kk == 100 )  EXIT
1433                ENDDO
1434             ENDIF
1435
1436             IF ( kk < 100  .AND.  uv_heights(kk+1) /= 9999999.9_wp )  THEN
1437                u_init(k) = u_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
1438                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
1439                                       ( u_profile(kk+1) - u_profile(kk) )
1440                v_init(k) = v_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
1441                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
1442                                       ( v_profile(kk+1) - v_profile(kk) )
1443             ELSE
1444                u_init(k) = u_profile(kk)
1445                v_init(k) = v_profile(kk)
1446             ENDIF
1447
1448          ENDDO
1449
1450       ELSE
1451
1452          message_string = 'u_profile(1) and v_profile(1) must be 0.0'
1453          CALL message( 'check_parameters', 'PA0346', 1, 2, 0, 6, 0 )
1454
1455       ENDIF
1456
1457!
1458!--    Compute initial temperature profile using the given temperature gradients
1459       IF (  .NOT.  neutral )  THEN
1460          CALL init_vertical_profiles( pt_vertical_gradient_level_ind,          &
1461                                       pt_vertical_gradient_level,              &
1462                                       pt_vertical_gradient, pt_init,           &
1463                                       pt_surface, bc_pt_t_val )
1464       ENDIF
1465!
1466!--    Compute initial humidity profile using the given humidity gradients
1467       IF ( humidity )  THEN
1468          CALL init_vertical_profiles( q_vertical_gradient_level_ind,          &
1469                                       q_vertical_gradient_level,              &
1470                                       q_vertical_gradient, q_init,            &
1471                                       q_surface, bc_q_t_val )
1472       ENDIF
1473!
1474!--    Compute initial scalar profile using the given scalar gradients
1475       IF ( passive_scalar )  THEN
1476          CALL init_vertical_profiles( s_vertical_gradient_level_ind,          &
1477                                       s_vertical_gradient_level,              &
1478                                       s_vertical_gradient, s_init,            &
1479                                       s_surface, bc_s_t_val )
1480       ENDIF
1481!
1482!--    If required, compute initial salinity profile using the given salinity
1483!--    gradients
1484       IF ( ocean )  THEN
1485          CALL init_vertical_profiles( sa_vertical_gradient_level_ind,          &
1486                                       sa_vertical_gradient_level,              &
1487                                       sa_vertical_gradient, sa_init,           &
1488                                       sa_surface, dum )
1489       ENDIF
1490
1491
1492    ENDIF
1493
1494!
1495!-- Check if the control parameter use_subsidence_tendencies is used correctly
1496    IF ( use_subsidence_tendencies  .AND.  .NOT.  large_scale_subsidence )  THEN
1497       message_string = 'The usage of use_subsidence_tendencies ' //           &
1498                            'requires large_scale_subsidence = .T..'
1499       CALL message( 'check_parameters', 'PA0396', 1, 2, 0, 6, 0 )
1500    ELSEIF ( use_subsidence_tendencies  .AND.  .NOT. large_scale_forcing )  THEN
1501       message_string = 'The usage of use_subsidence_tendencies ' //           &
1502                            'requires large_scale_forcing = .T..'
1503       CALL message( 'check_parameters', 'PA0397', 1, 2, 0, 6, 0 )
1504    ENDIF
1505
1506!
1507!-- Initialize large scale subsidence if required
1508    If ( large_scale_subsidence )  THEN
1509       IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp  .AND.            &
1510                                     .NOT.  large_scale_forcing )  THEN
1511          CALL init_w_subsidence
1512       ENDIF
1513!
1514!--    In case large_scale_forcing is used, profiles for subsidence velocity
1515!--    are read in from file LSF_DATA
1516
1517       IF ( subs_vertical_gradient_level(1) == -9999999.9_wp  .AND.            &
1518            .NOT.  large_scale_forcing )  THEN
1519          message_string = 'There is no default large scale vertical ' //      &
1520                           'velocity profile set. Specify the subsidence ' //  &
1521                           'velocity profile via subs_vertical_gradient ' //   &
1522                           'and subs_vertical_gradient_level.'
1523          CALL message( 'check_parameters', 'PA0380', 1, 2, 0, 6, 0 )
1524       ENDIF
1525    ELSE
1526        IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp )  THEN
1527           message_string = 'Enable usage of large scale subsidence by ' //    &
1528                            'setting large_scale_subsidence = .T..'
1529          CALL message( 'check_parameters', 'PA0381', 1, 2, 0, 6, 0 )
1530        ENDIF
1531    ENDIF
1532
1533!
1534!-- Compute Coriolis parameter
1535    f  = 2.0_wp * omega * SIN( phi / 180.0_wp * pi )
1536    fs = 2.0_wp * omega * COS( phi / 180.0_wp * pi )
1537
1538!
1539!-- Check and set buoyancy related parameters and switches
1540    IF ( reference_state == 'horizontal_average' )  THEN
1541       CONTINUE
1542    ELSEIF ( reference_state == 'initial_profile' )  THEN
1543       use_initial_profile_as_reference = .TRUE.
1544    ELSEIF ( reference_state == 'single_value' )  THEN
1545       use_single_reference_value = .TRUE.
1546       IF ( pt_reference == 9999999.9_wp )  pt_reference = pt_surface
1547       vpt_reference = pt_reference * ( 1.0_wp + 0.61_wp * q_surface )
1548    ELSE
1549       message_string = 'illegal value for reference_state: "' //              &
1550                        TRIM( reference_state ) // '"'
1551       CALL message( 'check_parameters', 'PA0056', 1, 2, 0, 6, 0 )
1552    ENDIF
1553
1554!
1555!-- Sign of buoyancy/stability terms
1556    IF ( ocean )  atmos_ocean_sign = -1.0_wp
1557
1558!
1559!-- Ocean version must use flux boundary conditions at the top
1560    IF ( ocean .AND. .NOT. use_top_fluxes )  THEN
1561       message_string = 'use_top_fluxes must be .TRUE. in ocean mode'
1562       CALL message( 'check_parameters', 'PA0042', 1, 2, 0, 6, 0 )
1563    ENDIF
1564
1565!
1566!-- In case of a given slope, compute the relevant quantities
1567    IF ( alpha_surface /= 0.0_wp )  THEN
1568       IF ( ABS( alpha_surface ) > 90.0_wp )  THEN
1569          WRITE( message_string, * ) 'ABS( alpha_surface = ', alpha_surface,   &
1570                                     ' ) must be < 90.0'
1571          CALL message( 'check_parameters', 'PA0043', 1, 2, 0, 6, 0 )
1572       ENDIF
1573       sloping_surface = .TRUE.
1574       cos_alpha_surface = COS( alpha_surface / 180.0_wp * pi )
1575       sin_alpha_surface = SIN( alpha_surface / 180.0_wp * pi )
1576    ENDIF
1577
1578!
1579!-- Check time step and cfl_factor
1580    IF ( dt /= -1.0_wp )  THEN
1581       IF ( dt <= 0.0_wp  .AND.  dt /= -1.0_wp )  THEN
1582          WRITE( message_string, * ) 'dt = ', dt , ' <= 0.0'
1583          CALL message( 'check_parameters', 'PA0044', 1, 2, 0, 6, 0 )
1584       ENDIF
1585       dt_3d = dt
1586       dt_fixed = .TRUE.
1587    ENDIF
1588
1589    IF ( cfl_factor <= 0.0_wp  .OR.  cfl_factor > 1.0_wp )  THEN
1590       IF ( cfl_factor == -1.0_wp )  THEN
1591          IF ( timestep_scheme == 'runge-kutta-2' )  THEN
1592             cfl_factor = 0.8_wp
1593          ELSEIF ( timestep_scheme == 'runge-kutta-3' )  THEN
1594             cfl_factor = 0.9_wp
1595          ELSE
1596             cfl_factor = 0.9_wp
1597          ENDIF
1598       ELSE
1599          WRITE( message_string, * ) 'cfl_factor = ', cfl_factor,              &
1600                 ' out of range & 0.0 < cfl_factor <= 1.0 is required'
1601          CALL message( 'check_parameters', 'PA0045', 1, 2, 0, 6, 0 )
1602       ENDIF
1603    ENDIF
1604
1605!
1606!-- Store simulated time at begin
1607    simulated_time_at_begin = simulated_time
1608
1609!
1610!-- Store reference time for coupled runs and change the coupling flag,
1611!-- if ...
1612    IF ( simulated_time == 0.0_wp )  THEN
1613       IF ( coupling_start_time == 0.0_wp )  THEN
1614          time_since_reference_point = 0.0_wp
1615       ELSEIF ( time_since_reference_point < 0.0_wp )  THEN
1616          run_coupled = .FALSE.
1617       ENDIF
1618    ENDIF
1619
1620!
1621!-- Set wind speed in the Galilei-transformed system
1622    IF ( galilei_transformation )  THEN
1623       IF ( use_ug_for_galilei_tr                    .AND.                     &
1624            ug_vertical_gradient_level(1) == 0.0_wp  .AND.                     &
1625            ug_vertical_gradient(1) == 0.0_wp        .AND.                     &
1626            vg_vertical_gradient_level(1) == 0.0_wp  .AND.                     &
1627            vg_vertical_gradient(1) == 0.0_wp )  THEN
1628          u_gtrans = ug_surface * 0.6_wp
1629          v_gtrans = vg_surface * 0.6_wp
1630       ELSEIF ( use_ug_for_galilei_tr  .AND.                                   &
1631                ( ug_vertical_gradient_level(1) /= 0.0_wp  .OR.                &
1632                ug_vertical_gradient(1) /= 0.0_wp ) )  THEN
1633          message_string = 'baroclinity (ug) not allowed simultaneously' //    &
1634                           ' with galilei transformation'
1635          CALL message( 'check_parameters', 'PA0046', 1, 2, 0, 6, 0 )
1636       ELSEIF ( use_ug_for_galilei_tr  .AND.                                   &
1637                ( vg_vertical_gradient_level(1) /= 0.0_wp  .OR.                &
1638                vg_vertical_gradient(1) /= 0.0_wp ) )  THEN
1639          message_string = 'baroclinity (vg) not allowed simultaneously' //    &
1640                           ' with galilei transformation'
1641          CALL message( 'check_parameters', 'PA0047', 1, 2, 0, 6, 0 )
1642       ELSE
1643          message_string = 'variable translation speed used for galilei-' //   &
1644             'transformation, which may cause & instabilities in stably ' //   &
1645             'stratified regions'
1646          CALL message( 'check_parameters', 'PA0048', 0, 1, 0, 6, 0 )
1647       ENDIF
1648    ENDIF
1649
1650!
1651!-- In case of using a prandtl-layer, calculated (or prescribed) surface
1652!-- fluxes have to be used in the diffusion-terms
1653    IF ( constant_flux_layer )  use_surface_fluxes = .TRUE.
1654
1655!
1656!-- Check boundary conditions and set internal variables:
1657!-- Attention: the lateral boundary conditions have been already checked in
1658!-- parin
1659!
1660!-- Non-cyclic lateral boundaries require the multigrid method and Piascek-
1661!-- Willimas or Wicker - Skamarock advection scheme. Several schemes
1662!-- and tools do not work with non-cyclic boundary conditions.
1663    IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
1664       IF ( psolver(1:9) /= 'multigrid' )  THEN
1665          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1666                           'psolver = "' // TRIM( psolver ) // '"'
1667          CALL message( 'check_parameters', 'PA0051', 1, 2, 0, 6, 0 )
1668       ENDIF
1669       IF ( momentum_advec /= 'pw-scheme'  .AND.                               &
1670            momentum_advec /= 'ws-scheme' )  THEN
1671
1672          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1673                           'momentum_advec = "' // TRIM( momentum_advec ) // '"'
1674          CALL message( 'check_parameters', 'PA0052', 1, 2, 0, 6, 0 )
1675       ENDIF
1676       IF ( scalar_advec /= 'pw-scheme'  .AND.                                 &
1677            scalar_advec /= 'ws-scheme' )  THEN
1678          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1679                           'scalar_advec = "' // TRIM( scalar_advec ) // '"'
1680          CALL message( 'check_parameters', 'PA0053', 1, 2, 0, 6, 0 )
1681       ENDIF
1682       IF ( galilei_transformation )  THEN
1683          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1684                           'galilei_transformation = .T.'
1685          CALL message( 'check_parameters', 'PA0054', 1, 2, 0, 6, 0 )
1686       ENDIF
1687    ENDIF
1688
1689!
1690!-- Bottom boundary condition for the turbulent Kinetic energy
1691    IF ( bc_e_b == 'neumann' )  THEN
1692       ibc_e_b = 1
1693    ELSEIF ( bc_e_b == '(u*)**2+neumann' )  THEN
1694       ibc_e_b = 2
1695       IF ( .NOT. constant_flux_layer )  THEN
1696          bc_e_b = 'neumann'
1697          ibc_e_b = 1
1698          message_string = 'boundary condition bc_e_b changed to "' //         &
1699                           TRIM( bc_e_b ) // '"'
1700          CALL message( 'check_parameters', 'PA0057', 0, 1, 0, 6, 0 )
1701       ENDIF
1702    ELSE
1703       message_string = 'unknown boundary condition: bc_e_b = "' //            &
1704                        TRIM( bc_e_b ) // '"'
1705       CALL message( 'check_parameters', 'PA0058', 1, 2, 0, 6, 0 )
1706    ENDIF
1707
1708!
1709!-- Boundary conditions for perturbation pressure
1710    IF ( bc_p_b == 'dirichlet' )  THEN
1711       ibc_p_b = 0
1712    ELSEIF ( bc_p_b == 'neumann' )  THEN
1713       ibc_p_b = 1
1714    ELSE
1715       message_string = 'unknown boundary condition: bc_p_b = "' //            &
1716                        TRIM( bc_p_b ) // '"'
1717       CALL message( 'check_parameters', 'PA0059', 1, 2, 0, 6, 0 )
1718    ENDIF
1719
1720    IF ( bc_p_t == 'dirichlet' )  THEN
1721       ibc_p_t = 0
1722!-- TO_DO: later set bc_p_t to neumann before, in case of nested domain
1723    ELSEIF ( bc_p_t == 'neumann' .OR. bc_p_t == 'nested' )  THEN
1724       ibc_p_t = 1
1725    ELSE
1726       message_string = 'unknown boundary condition: bc_p_t = "' //            &
1727                        TRIM( bc_p_t ) // '"'
1728       CALL message( 'check_parameters', 'PA0061', 1, 2, 0, 6, 0 )
1729    ENDIF
1730
1731!
1732!-- Boundary conditions for potential temperature
1733    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
1734       ibc_pt_b = 2
1735    ELSE
1736       IF ( bc_pt_b == 'dirichlet' )  THEN
1737          ibc_pt_b = 0
1738       ELSEIF ( bc_pt_b == 'neumann' )  THEN
1739          ibc_pt_b = 1
1740       ELSE
1741          message_string = 'unknown boundary condition: bc_pt_b = "' //        &
1742                           TRIM( bc_pt_b ) // '"'
1743          CALL message( 'check_parameters', 'PA0062', 1, 2, 0, 6, 0 )
1744       ENDIF
1745    ENDIF
1746
1747    IF ( bc_pt_t == 'dirichlet' )  THEN
1748       ibc_pt_t = 0
1749    ELSEIF ( bc_pt_t == 'neumann' )  THEN
1750       ibc_pt_t = 1
1751    ELSEIF ( bc_pt_t == 'initial_gradient' )  THEN
1752       ibc_pt_t = 2
1753    ELSEIF ( bc_pt_t == 'nested' )  THEN
1754       ibc_pt_t = 3
1755    ELSE
1756       message_string = 'unknown boundary condition: bc_pt_t = "' //           &
1757                        TRIM( bc_pt_t ) // '"'
1758       CALL message( 'check_parameters', 'PA0063', 1, 2, 0, 6, 0 )
1759    ENDIF
1760
1761    IF ( ANY( wall_heatflux /= 0.0_wp )  .AND.                        &
1762         surface_heatflux == 9999999.9_wp )  THEN
1763       message_string = 'wall_heatflux additionally requires ' //     &
1764                        'setting of surface_heatflux'
1765       CALL message( 'check_parameters', 'PA0443', 1, 2, 0, 6, 0 )
1766    ENDIF
1767
1768!
1769!   This IF clause needs revision, got too complex!!
1770    IF ( surface_heatflux == 9999999.9_wp  )  THEN
1771       constant_heatflux = .FALSE.
1772       IF ( large_scale_forcing  .OR.  land_surface  .OR.  urban_surface )  THEN
1773          IF ( ibc_pt_b == 0 )  THEN
1774             constant_heatflux = .FALSE.
1775          ELSEIF ( ibc_pt_b == 1 )  THEN
1776             constant_heatflux = .TRUE.
1777             IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.   &
1778                  .NOT.  land_surface )  THEN
1779                surface_heatflux = shf_surf(1)
1780             ELSE
1781                surface_heatflux = 0.0_wp
1782             ENDIF
1783          ENDIF
1784       ENDIF
1785    ELSE
1786        constant_heatflux = .TRUE.
1787        IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.        &
1788             large_scale_forcing  .AND.  .NOT.  land_surface )  THEN
1789           surface_heatflux = shf_surf(1)
1790        ENDIF
1791    ENDIF
1792
1793    IF ( top_heatflux     == 9999999.9_wp )  constant_top_heatflux = .FALSE.
1794
1795    IF ( neutral )  THEN
1796
1797       IF ( surface_heatflux /= 0.0_wp  .AND.                                  &
1798            surface_heatflux /= 9999999.9_wp )  THEN
1799          message_string = 'heatflux must not be set for pure neutral flow'
1800          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
1801       ENDIF
1802
1803       IF ( top_heatflux /= 0.0_wp  .AND.  top_heatflux /= 9999999.9_wp )      &
1804       THEN
1805          message_string = 'heatflux must not be set for pure neutral flow'
1806          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
1807       ENDIF
1808
1809    ENDIF
1810
1811    IF ( top_momentumflux_u /= 9999999.9_wp  .AND.                             &
1812         top_momentumflux_v /= 9999999.9_wp )  THEN
1813       constant_top_momentumflux = .TRUE.
1814    ELSEIF (  .NOT. ( top_momentumflux_u == 9999999.9_wp  .AND.                &
1815           top_momentumflux_v == 9999999.9_wp ) )  THEN
1816       message_string = 'both, top_momentumflux_u AND top_momentumflux_v ' //  &
1817                        'must be set'
1818       CALL message( 'check_parameters', 'PA0064', 1, 2, 0, 6, 0 )
1819    ENDIF
1820
1821!
1822!-- A given surface temperature implies Dirichlet boundary condition for
1823!-- temperature. In this case specification of a constant heat flux is
1824!-- forbidden.
1825    IF ( ibc_pt_b == 0  .AND.  constant_heatflux  .AND.                        &
1826         surface_heatflux /= 0.0_wp )  THEN
1827       message_string = 'boundary_condition: bc_pt_b = "' // TRIM( bc_pt_b ) //&
1828                        '& is not allowed with constant_heatflux = .TRUE.'
1829       CALL message( 'check_parameters', 'PA0065', 1, 2, 0, 6, 0 )
1830    ENDIF
1831    IF ( constant_heatflux  .AND.  pt_surface_initial_change /= 0.0_wp )  THEN
1832       WRITE ( message_string, * )  'constant_heatflux = .TRUE. is not allo',  &
1833               'wed with pt_surface_initial_change (/=0) = ',                  &
1834               pt_surface_initial_change
1835       CALL message( 'check_parameters', 'PA0066', 1, 2, 0, 6, 0 )
1836    ENDIF
1837
1838!
1839!-- A given temperature at the top implies Dirichlet boundary condition for
1840!-- temperature. In this case specification of a constant heat flux is
1841!-- forbidden.
1842    IF ( ibc_pt_t == 0  .AND.  constant_top_heatflux  .AND.                    &
1843         top_heatflux /= 0.0_wp )  THEN
1844       message_string = 'boundary_condition: bc_pt_t = "' // TRIM( bc_pt_t ) //&
1845                        '" is not allowed with constant_top_heatflux = .TRUE.'
1846       CALL message( 'check_parameters', 'PA0067', 1, 2, 0, 6, 0 )
1847    ENDIF
1848
1849!
1850!-- Boundary conditions for salinity
1851    IF ( ocean )  THEN
1852       IF ( bc_sa_t == 'dirichlet' )  THEN
1853          ibc_sa_t = 0
1854       ELSEIF ( bc_sa_t == 'neumann' )  THEN
1855          ibc_sa_t = 1
1856       ELSE
1857          message_string = 'unknown boundary condition: bc_sa_t = "' //        &
1858                           TRIM( bc_sa_t ) // '"'
1859          CALL message( 'check_parameters', 'PA0068', 1, 2, 0, 6, 0 )
1860       ENDIF
1861
1862       IF ( top_salinityflux == 9999999.9_wp )  constant_top_salinityflux = .FALSE.
1863       IF ( ibc_sa_t == 1  .AND.  top_salinityflux == 9999999.9_wp )  THEN
1864          message_string = 'boundary condition: bc_sa_t = "' //                &
1865                           TRIM( bc_sa_t ) // '" requires to set ' //          &
1866                           'top_salinityflux'
1867          CALL message( 'check_parameters', 'PA0069', 1, 2, 0, 6, 0 )
1868       ENDIF
1869
1870!
1871!--    A fixed salinity at the top implies Dirichlet boundary condition for
1872!--    salinity. In this case specification of a constant salinity flux is
1873!--    forbidden.
1874       IF ( ibc_sa_t == 0  .AND.  constant_top_salinityflux  .AND.             &
1875            top_salinityflux /= 0.0_wp )  THEN
1876          message_string = 'boundary condition: bc_sa_t = "' //                &
1877                           TRIM( bc_sa_t ) // '" is not allowed with ' //      &
1878                           'top_salinityflux /= 0.0'
1879          CALL message( 'check_parameters', 'PA0070', 1, 2, 0, 6, 0 )
1880       ENDIF
1881
1882    ENDIF
1883
1884!
1885!-- Set boundary conditions for total water content
1886    IF ( humidity )  THEN
1887
1888       IF ( ANY( wall_humidityflux /= 0.0_wp )  .AND.                        &
1889            surface_waterflux == 9999999.9_wp )  THEN
1890          message_string = 'wall_humidityflux additionally requires ' //     &
1891                           'setting of surface_waterflux'
1892          CALL message( 'check_parameters', 'PA0444', 1, 2, 0, 6, 0 )
1893       ENDIF
1894
1895       CALL set_bc_scalars( 'q', bc_q_b, bc_q_t, ibc_q_b, ibc_q_t,           &
1896                            'PA0071', 'PA0072' )
1897
1898       IF ( surface_waterflux == 9999999.9_wp  )  THEN
1899          constant_waterflux = .FALSE.
1900          IF ( large_scale_forcing .OR. land_surface )  THEN
1901             IF ( ibc_q_b == 0 )  THEN
1902                constant_waterflux = .FALSE.
1903             ELSEIF ( ibc_q_b == 1 )  THEN
1904                constant_waterflux = .TRUE.
1905                IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1906                   surface_waterflux = qsws_surf(1)
1907                ENDIF
1908             ENDIF
1909          ENDIF
1910       ELSE
1911          constant_waterflux = .TRUE.
1912          IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.      &
1913               large_scale_forcing )  THEN
1914             surface_waterflux = qsws_surf(1)
1915          ENDIF
1916       ENDIF
1917
1918       CALL check_bc_scalars( 'q', bc_q_b, ibc_q_b, 'PA0073', 'PA0074',        &
1919                              constant_waterflux, q_surface_initial_change )
1920
1921    ENDIF
1922
1923    IF ( passive_scalar )  THEN
1924
1925       IF ( ANY( wall_scalarflux /= 0.0_wp )  .AND.                        &
1926            surface_scalarflux == 9999999.9_wp )  THEN
1927          message_string = 'wall_scalarflux additionally requires ' //     &
1928                           'setting of surface_scalarflux'
1929          CALL message( 'check_parameters', 'PA0445', 1, 2, 0, 6, 0 )
1930       ENDIF
1931
1932       IF ( surface_scalarflux == 9999999.9_wp )  constant_scalarflux = .FALSE.
1933
1934       CALL set_bc_scalars( 's', bc_s_b, bc_s_t, ibc_s_b, ibc_s_t,             &
1935                            'PA0071', 'PA0072' )
1936
1937       CALL check_bc_scalars( 's', bc_s_b, ibc_s_b, 'PA0073', 'PA0074',        &
1938                              constant_scalarflux, s_surface_initial_change )
1939
1940       IF ( top_scalarflux == 9999999.9_wp )  constant_top_scalarflux = .FALSE.
1941!
1942!--    A fixed scalar concentration at the top implies Dirichlet boundary
1943!--    condition for scalar. Hence, in this case specification of a constant
1944!--    scalar flux is forbidden.
1945       IF ( ( ibc_s_t == 0 .OR. ibc_s_t == 2 )  .AND.  constant_top_scalarflux &
1946               .AND.  top_scalarflux /= 0.0_wp )  THEN
1947          message_string = 'boundary condition: bc_s_t = "' //                 &
1948                           TRIM( bc_sa_t ) // '" is not allowed with ' //      &
1949                           'top_scalarflux /= 0.0'
1950          CALL message( 'check_parameters', 'PA0441', 1, 2, 0, 6, 0 )
1951       ENDIF
1952    ENDIF
1953!
1954!-- Boundary conditions for horizontal components of wind speed
1955    IF ( bc_uv_b == 'dirichlet' )  THEN
1956       ibc_uv_b = 0
1957    ELSEIF ( bc_uv_b == 'neumann' )  THEN
1958       ibc_uv_b = 1
1959       IF ( constant_flux_layer )  THEN
1960          message_string = 'boundary condition: bc_uv_b = "' //                &
1961               TRIM( bc_uv_b ) // '" is not allowed with constant_flux_layer'  &
1962               // ' = .TRUE.'
1963          CALL message( 'check_parameters', 'PA0075', 1, 2, 0, 6, 0 )
1964       ENDIF
1965    ELSE
1966       message_string = 'unknown boundary condition: bc_uv_b = "' //           &
1967                        TRIM( bc_uv_b ) // '"'
1968       CALL message( 'check_parameters', 'PA0076', 1, 2, 0, 6, 0 )
1969    ENDIF
1970!
1971!-- In case of coupled simulations u and v at the ground in atmosphere will be
1972!-- assigned with the u and v values of the ocean surface
1973    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
1974       ibc_uv_b = 2
1975    ENDIF
1976
1977    IF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
1978       bc_uv_t = 'neumann'
1979       ibc_uv_t = 1
1980    ELSE
1981       IF ( bc_uv_t == 'dirichlet' .OR. bc_uv_t == 'dirichlet_0' )  THEN
1982          ibc_uv_t = 0
1983          IF ( bc_uv_t == 'dirichlet_0' )  THEN
1984!
1985!--          Velocities for the initial u,v-profiles are set zero at the top
1986!--          in case of dirichlet_0 conditions
1987             u_init(nzt+1)    = 0.0_wp
1988             v_init(nzt+1)    = 0.0_wp
1989          ENDIF
1990       ELSEIF ( bc_uv_t == 'neumann' )  THEN
1991          ibc_uv_t = 1
1992       ELSEIF ( bc_uv_t == 'nested' )  THEN
1993          ibc_uv_t = 3
1994       ELSE
1995          message_string = 'unknown boundary condition: bc_uv_t = "' //        &
1996                           TRIM( bc_uv_t ) // '"'
1997          CALL message( 'check_parameters', 'PA0077', 1, 2, 0, 6, 0 )
1998       ENDIF
1999    ENDIF
2000
2001!
2002!-- Compute and check, respectively, the Rayleigh Damping parameter
2003    IF ( rayleigh_damping_factor == -1.0_wp )  THEN
2004       rayleigh_damping_factor = 0.0_wp
2005    ELSE
2006       IF ( rayleigh_damping_factor < 0.0_wp  .OR.                             &
2007            rayleigh_damping_factor > 1.0_wp )  THEN
2008          WRITE( message_string, * )  'rayleigh_damping_factor = ',            &
2009                              rayleigh_damping_factor, ' out of range [0.0,1.0]'
2010          CALL message( 'check_parameters', 'PA0078', 1, 2, 0, 6, 0 )
2011       ENDIF
2012    ENDIF
2013
2014    IF ( rayleigh_damping_height == -1.0_wp )  THEN
2015       IF (  .NOT.  ocean )  THEN
2016          rayleigh_damping_height = 0.66666666666_wp * zu(nzt)
2017       ELSE
2018          rayleigh_damping_height = 0.66666666666_wp * zu(nzb)
2019       ENDIF
2020    ELSE
2021       IF (  .NOT.  ocean )  THEN
2022          IF ( rayleigh_damping_height < 0.0_wp  .OR.                          &
2023               rayleigh_damping_height > zu(nzt) )  THEN
2024             WRITE( message_string, * )  'rayleigh_damping_height = ',         &
2025                   rayleigh_damping_height, ' out of range [0.0,', zu(nzt), ']'
2026             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
2027          ENDIF
2028       ELSE
2029          IF ( rayleigh_damping_height > 0.0_wp  .OR.                          &
2030               rayleigh_damping_height < zu(nzb) )  THEN
2031             WRITE( message_string, * )  'rayleigh_damping_height = ',         &
2032                   rayleigh_damping_height, ' out of range [0.0,', zu(nzb), ']'
2033             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
2034          ENDIF
2035       ENDIF
2036    ENDIF
2037
2038!
2039!-- Check number of chosen statistic regions
2040    IF ( statistic_regions < 0 )  THEN
2041       WRITE ( message_string, * ) 'number of statistic_regions = ',           &
2042                   statistic_regions+1, ' is not allowed'
2043       CALL message( 'check_parameters', 'PA0082', 1, 2, 0, 6, 0 )
2044    ENDIF
2045    IF ( normalizing_region > statistic_regions  .OR.                          &
2046         normalizing_region < 0)  THEN
2047       WRITE ( message_string, * ) 'normalizing_region = ',                    &
2048                normalizing_region, ' must be >= 0 and <= ',statistic_regions, &
2049                ' (value of statistic_regions)'
2050       CALL message( 'check_parameters', 'PA0083', 1, 2, 0, 6, 0 )
2051    ENDIF
2052
2053!
2054!-- Set the default intervals for data output, if necessary
2055!-- NOTE: dt_dosp has already been set in spectra_parin
2056    IF ( dt_data_output /= 9999999.9_wp )  THEN
2057       IF ( dt_dopr           == 9999999.9_wp )  dt_dopr           = dt_data_output
2058       IF ( dt_dopts          == 9999999.9_wp )  dt_dopts          = dt_data_output
2059       IF ( dt_do2d_xy        == 9999999.9_wp )  dt_do2d_xy        = dt_data_output
2060       IF ( dt_do2d_xz        == 9999999.9_wp )  dt_do2d_xz        = dt_data_output
2061       IF ( dt_do2d_yz        == 9999999.9_wp )  dt_do2d_yz        = dt_data_output
2062       IF ( dt_do3d           == 9999999.9_wp )  dt_do3d           = dt_data_output
2063       IF ( dt_data_output_av == 9999999.9_wp )  dt_data_output_av = dt_data_output
2064       DO  mid = 1, max_masks
2065          IF ( dt_domask(mid) == 9999999.9_wp )  dt_domask(mid)    = dt_data_output
2066       ENDDO
2067    ENDIF
2068
2069!
2070!-- Set the default skip time intervals for data output, if necessary
2071    IF ( skip_time_dopr    == 9999999.9_wp )                                   &
2072                                       skip_time_dopr    = skip_time_data_output
2073    IF ( skip_time_do2d_xy == 9999999.9_wp )                                   &
2074                                       skip_time_do2d_xy = skip_time_data_output
2075    IF ( skip_time_do2d_xz == 9999999.9_wp )                                   &
2076                                       skip_time_do2d_xz = skip_time_data_output
2077    IF ( skip_time_do2d_yz == 9999999.9_wp )                                   &
2078                                       skip_time_do2d_yz = skip_time_data_output
2079    IF ( skip_time_do3d    == 9999999.9_wp )                                   &
2080                                       skip_time_do3d    = skip_time_data_output
2081    IF ( skip_time_data_output_av == 9999999.9_wp )                            &
2082                                skip_time_data_output_av = skip_time_data_output
2083    DO  mid = 1, max_masks
2084       IF ( skip_time_domask(mid) == 9999999.9_wp )                            &
2085                                skip_time_domask(mid)    = skip_time_data_output
2086    ENDDO
2087
2088!
2089!-- Check the average intervals (first for 3d-data, then for profiles)
2090    IF ( averaging_interval > dt_data_output_av )  THEN
2091       WRITE( message_string, * )  'averaging_interval = ',                    &
2092             averaging_interval, ' must be <= dt_data_output = ', dt_data_output
2093       CALL message( 'check_parameters', 'PA0085', 1, 2, 0, 6, 0 )
2094    ENDIF
2095
2096    IF ( averaging_interval_pr == 9999999.9_wp )  THEN
2097       averaging_interval_pr = averaging_interval
2098    ENDIF
2099
2100    IF ( averaging_interval_pr > dt_dopr )  THEN
2101       WRITE( message_string, * )  'averaging_interval_pr = ',                 &
2102             averaging_interval_pr, ' must be <= dt_dopr = ', dt_dopr
2103       CALL message( 'check_parameters', 'PA0086', 1, 2, 0, 6, 0 )
2104    ENDIF
2105
2106!
2107!-- Set the default interval for profiles entering the temporal average
2108    IF ( dt_averaging_input_pr == 9999999.9_wp )  THEN
2109       dt_averaging_input_pr = dt_averaging_input
2110    ENDIF
2111
2112!
2113!-- Set the default interval for the output of timeseries to a reasonable
2114!-- value (tries to minimize the number of calls of flow_statistics)
2115    IF ( dt_dots == 9999999.9_wp )  THEN
2116       IF ( averaging_interval_pr == 0.0_wp )  THEN
2117          dt_dots = MIN( dt_run_control, dt_dopr )
2118       ELSE
2119          dt_dots = MIN( dt_run_control, dt_averaging_input_pr )
2120       ENDIF
2121    ENDIF
2122
2123!
2124!-- Check the sample rate for averaging (first for 3d-data, then for profiles)
2125    IF ( dt_averaging_input > averaging_interval )  THEN
2126       WRITE( message_string, * )  'dt_averaging_input = ',                    &
2127                dt_averaging_input, ' must be <= averaging_interval = ',       &
2128                averaging_interval
2129       CALL message( 'check_parameters', 'PA0088', 1, 2, 0, 6, 0 )
2130    ENDIF
2131
2132    IF ( dt_averaging_input_pr > averaging_interval_pr )  THEN
2133       WRITE( message_string, * )  'dt_averaging_input_pr = ',                 &
2134                dt_averaging_input_pr, ' must be <= averaging_interval_pr = ', &
2135                averaging_interval_pr
2136       CALL message( 'check_parameters', 'PA0089', 1, 2, 0, 6, 0 )
2137    ENDIF
2138
2139!
2140!-- Set the default value for the integration interval of precipitation amount
2141    IF ( microphysics_seifert  .OR.  microphysics_kessler )  THEN
2142       IF ( precipitation_amount_interval == 9999999.9_wp )  THEN
2143          precipitation_amount_interval = dt_do2d_xy
2144       ELSE
2145          IF ( precipitation_amount_interval > dt_do2d_xy )  THEN
2146             WRITE( message_string, * )  'precipitation_amount_interval = ',   &
2147                 precipitation_amount_interval, ' must not be larger than ',   &
2148                 'dt_do2d_xy = ', dt_do2d_xy
2149             CALL message( 'check_parameters', 'PA0090', 1, 2, 0, 6, 0 )
2150          ENDIF
2151       ENDIF
2152    ENDIF
2153
2154!
2155!-- Determine the number of output profiles and check whether they are
2156!-- permissible
2157    DO  WHILE ( data_output_pr(dopr_n+1) /= '          ' )
2158
2159       dopr_n = dopr_n + 1
2160       i = dopr_n
2161
2162!
2163!--    Determine internal profile number (for hom, homs)
2164!--    and store height levels
2165       SELECT CASE ( TRIM( data_output_pr(i) ) )
2166
2167          CASE ( 'u', '#u' )
2168             dopr_index(i) = 1
2169             dopr_unit(i)  = 'm/s'
2170             hom(:,2,1,:)  = SPREAD( zu, 2, statistic_regions+1 )
2171             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2172                dopr_initial_index(i) = 5
2173                hom(:,2,5,:)          = SPREAD( zu, 2, statistic_regions+1 )
2174                data_output_pr(i)     = data_output_pr(i)(2:)
2175             ENDIF
2176
2177          CASE ( 'v', '#v' )
2178             dopr_index(i) = 2
2179             dopr_unit(i)  = 'm/s'
2180             hom(:,2,2,:)  = SPREAD( zu, 2, statistic_regions+1 )
2181             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2182                dopr_initial_index(i) = 6
2183                hom(:,2,6,:)          = SPREAD( zu, 2, statistic_regions+1 )
2184                data_output_pr(i)     = data_output_pr(i)(2:)
2185             ENDIF
2186
2187          CASE ( 'w' )
2188             dopr_index(i) = 3
2189             dopr_unit(i)  = 'm/s'
2190             hom(:,2,3,:)  = SPREAD( zw, 2, statistic_regions+1 )
2191
2192          CASE ( 'pt', '#pt' )
2193             IF ( .NOT. cloud_physics ) THEN
2194                dopr_index(i) = 4
2195                dopr_unit(i)  = 'K'
2196                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
2197                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2198                   dopr_initial_index(i) = 7
2199                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
2200                   hom(nzb,2,7,:)        = 0.0_wp    ! because zu(nzb) is negative
2201                   data_output_pr(i)     = data_output_pr(i)(2:)
2202                ENDIF
2203             ELSE
2204                dopr_index(i) = 43
2205                dopr_unit(i)  = 'K'
2206                hom(:,2,43,:)  = SPREAD( zu, 2, statistic_regions+1 )
2207                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2208                   dopr_initial_index(i) = 28
2209                   hom(:,2,28,:)         = SPREAD( zu, 2, statistic_regions+1 )
2210                   hom(nzb,2,28,:)       = 0.0_wp    ! because zu(nzb) is negative
2211                   data_output_pr(i)     = data_output_pr(i)(2:)
2212                ENDIF
2213             ENDIF
2214
2215          CASE ( 'e' )
2216             dopr_index(i)  = 8
2217             dopr_unit(i)   = 'm2/s2'
2218             hom(:,2,8,:)   = SPREAD( zu, 2, statistic_regions+1 )
2219             hom(nzb,2,8,:) = 0.0_wp
2220
2221          CASE ( 'km', '#km' )
2222             dopr_index(i)  = 9
2223             dopr_unit(i)   = 'm2/s'
2224             hom(:,2,9,:)   = SPREAD( zu, 2, statistic_regions+1 )
2225             hom(nzb,2,9,:) = 0.0_wp
2226             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2227                dopr_initial_index(i) = 23
2228                hom(:,2,23,:)         = hom(:,2,9,:)
2229                data_output_pr(i)     = data_output_pr(i)(2:)
2230             ENDIF
2231
2232          CASE ( 'kh', '#kh' )
2233             dopr_index(i)   = 10
2234             dopr_unit(i)    = 'm2/s'
2235             hom(:,2,10,:)   = SPREAD( zu, 2, statistic_regions+1 )
2236             hom(nzb,2,10,:) = 0.0_wp
2237             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2238                dopr_initial_index(i) = 24
2239                hom(:,2,24,:)         = hom(:,2,10,:)
2240                data_output_pr(i)     = data_output_pr(i)(2:)
2241             ENDIF
2242
2243          CASE ( 'l', '#l' )
2244             dopr_index(i)   = 11
2245             dopr_unit(i)    = 'm'
2246             hom(:,2,11,:)   = SPREAD( zu, 2, statistic_regions+1 )
2247             hom(nzb,2,11,:) = 0.0_wp
2248             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2249                dopr_initial_index(i) = 25
2250                hom(:,2,25,:)         = hom(:,2,11,:)
2251                data_output_pr(i)     = data_output_pr(i)(2:)
2252             ENDIF
2253
2254          CASE ( 'w"u"' )
2255             dopr_index(i) = 12
2256             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2257             hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 )
2258             IF ( constant_flux_layer )  hom(nzb,2,12,:) = zu(1)
2259
2260          CASE ( 'w*u*' )
2261             dopr_index(i) = 13
2262             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2263             hom(:,2,13,:) = SPREAD( zw, 2, statistic_regions+1 )
2264
2265          CASE ( 'w"v"' )
2266             dopr_index(i) = 14
2267             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2268             hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 )
2269             IF ( constant_flux_layer )  hom(nzb,2,14,:) = zu(1)
2270
2271          CASE ( 'w*v*' )
2272             dopr_index(i) = 15
2273             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2274             hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 )
2275
2276          CASE ( 'w"pt"' )
2277             dopr_index(i) = 16
2278             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2279             hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 )
2280
2281          CASE ( 'w*pt*' )
2282             dopr_index(i) = 17
2283             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2284             hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 )
2285
2286          CASE ( 'wpt' )
2287             dopr_index(i) = 18
2288             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2289             hom(:,2,18,:) = SPREAD( zw, 2, statistic_regions+1 )
2290
2291          CASE ( 'wu' )
2292             dopr_index(i) = 19
2293             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2294             hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 )
2295             IF ( constant_flux_layer )  hom(nzb,2,19,:) = zu(1)
2296
2297          CASE ( 'wv' )
2298             dopr_index(i) = 20
2299             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2300             hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 )
2301             IF ( constant_flux_layer )  hom(nzb,2,20,:) = zu(1)
2302
2303          CASE ( 'w*pt*BC' )
2304             dopr_index(i) = 21
2305             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2306             hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 )
2307
2308          CASE ( 'wptBC' )
2309             dopr_index(i) = 22
2310             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2311             hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
2312
2313          CASE ( 'sa', '#sa' )
2314             IF ( .NOT. ocean )  THEN
2315                message_string = 'data_output_pr = ' // &
2316                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2317                                 'lemented for ocean = .FALSE.'
2318                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2319             ELSE
2320                dopr_index(i) = 23
2321                dopr_unit(i)  = 'psu'
2322                hom(:,2,23,:) = SPREAD( zu, 2, statistic_regions+1 )
2323                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2324                   dopr_initial_index(i) = 26
2325                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2326                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2327                   data_output_pr(i)     = data_output_pr(i)(2:)
2328                ENDIF
2329             ENDIF
2330
2331          CASE ( 'u*2' )
2332             dopr_index(i) = 30
2333             dopr_unit(i)  = 'm2/s2'
2334             hom(:,2,30,:) = SPREAD( zu, 2, statistic_regions+1 )
2335
2336          CASE ( 'v*2' )
2337             dopr_index(i) = 31
2338             dopr_unit(i)  = 'm2/s2'
2339             hom(:,2,31,:) = SPREAD( zu, 2, statistic_regions+1 )
2340
2341          CASE ( 'w*2' )
2342             dopr_index(i) = 32
2343             dopr_unit(i)  = 'm2/s2'
2344             hom(:,2,32,:) = SPREAD( zw, 2, statistic_regions+1 )
2345
2346          CASE ( 'pt*2' )
2347             dopr_index(i) = 33
2348             dopr_unit(i)  = 'K2'
2349             hom(:,2,33,:) = SPREAD( zu, 2, statistic_regions+1 )
2350
2351          CASE ( 'e*' )
2352             dopr_index(i) = 34
2353             dopr_unit(i)  = 'm2/s2'
2354             hom(:,2,34,:) = SPREAD( zu, 2, statistic_regions+1 )
2355
2356          CASE ( 'w*2pt*' )
2357             dopr_index(i) = 35
2358             dopr_unit(i)  = 'K m2/s2'
2359             hom(:,2,35,:) = SPREAD( zw, 2, statistic_regions+1 )
2360
2361          CASE ( 'w*pt*2' )
2362             dopr_index(i) = 36
2363             dopr_unit(i)  = 'K2 m/s'
2364             hom(:,2,36,:) = SPREAD( zw, 2, statistic_regions+1 )
2365
2366          CASE ( 'w*e*' )
2367             dopr_index(i) = 37
2368             dopr_unit(i)  = 'm3/s3'
2369             hom(:,2,37,:) = SPREAD( zw, 2, statistic_regions+1 )
2370
2371          CASE ( 'w*3' )
2372             dopr_index(i) = 38
2373             dopr_unit(i)  = 'm3/s3'
2374             hom(:,2,38,:) = SPREAD( zw, 2, statistic_regions+1 )
2375
2376          CASE ( 'Sw' )
2377             dopr_index(i) = 39
2378             dopr_unit(i)  = 'none'
2379             hom(:,2,39,:) = SPREAD( zw, 2, statistic_regions+1 )
2380
2381          CASE ( 'p' )
2382             dopr_index(i) = 40
2383             dopr_unit(i)  = 'Pa'
2384             hom(:,2,40,:) = SPREAD( zu, 2, statistic_regions+1 )
2385
2386          CASE ( 'q', '#q' )
2387             IF ( .NOT. humidity )  THEN
2388                message_string = 'data_output_pr = ' // &
2389                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2390                                 'lemented for humidity = .FALSE.'
2391                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2392             ELSE
2393                dopr_index(i) = 41
2394                dopr_unit(i)  = 'kg/kg'
2395                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2396                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2397                   dopr_initial_index(i) = 26
2398                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2399                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2400                   data_output_pr(i)     = data_output_pr(i)(2:)
2401                ENDIF
2402             ENDIF
2403
2404          CASE ( 's', '#s' )
2405             IF ( .NOT. passive_scalar )  THEN
2406                message_string = 'data_output_pr = ' //                        &
2407                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2408                                 'lemented for passive_scalar = .FALSE.'
2409                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2410             ELSE
2411                dopr_index(i) = 117
2412                dopr_unit(i)  = 'kg/m3'
2413                hom(:,2,117,:) = SPREAD( zu, 2, statistic_regions+1 )
2414                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2415                   dopr_initial_index(i) = 117
2416                   hom(:,2,117,:)        = SPREAD( zu, 2, statistic_regions+1 )
2417                   hom(nzb,2,117,:)      = 0.0_wp    ! because zu(nzb) is negative
2418                   data_output_pr(i)     = data_output_pr(i)(2:)
2419                ENDIF
2420             ENDIF
2421
2422          CASE ( 'qv', '#qv' )
2423             IF ( .NOT. cloud_physics ) THEN
2424                dopr_index(i) = 41
2425                dopr_unit(i)  = 'kg/kg'
2426                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2427                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2428                   dopr_initial_index(i) = 26
2429                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2430                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2431                   data_output_pr(i)     = data_output_pr(i)(2:)
2432                ENDIF
2433             ELSE
2434                dopr_index(i) = 42
2435                dopr_unit(i)  = 'kg/kg'
2436                hom(:,2,42,:) = SPREAD( zu, 2, statistic_regions+1 )
2437                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2438                   dopr_initial_index(i) = 27
2439                   hom(:,2,27,:)         = SPREAD( zu, 2, statistic_regions+1 )
2440                   hom(nzb,2,27,:)       = 0.0_wp   ! because zu(nzb) is negative
2441                   data_output_pr(i)     = data_output_pr(i)(2:)
2442                ENDIF
2443             ENDIF
2444
2445          CASE ( 'lpt', '#lpt' )
2446             IF ( .NOT. cloud_physics ) THEN
2447                message_string = 'data_output_pr = ' //                        &
2448                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2449                                 'lemented for cloud_physics = .FALSE.'
2450                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2451             ELSE
2452                dopr_index(i) = 4
2453                dopr_unit(i)  = 'K'
2454                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
2455                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2456                   dopr_initial_index(i) = 7
2457                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
2458                   hom(nzb,2,7,:)        = 0.0_wp    ! because zu(nzb) is negative
2459                   data_output_pr(i)     = data_output_pr(i)(2:)
2460                ENDIF
2461             ENDIF
2462
2463          CASE ( 'vpt', '#vpt' )
2464             dopr_index(i) = 44
2465             dopr_unit(i)  = 'K'
2466             hom(:,2,44,:) = SPREAD( zu, 2, statistic_regions+1 )
2467             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2468                dopr_initial_index(i) = 29
2469                hom(:,2,29,:)         = SPREAD( zu, 2, statistic_regions+1 )
2470                hom(nzb,2,29,:)       = 0.0_wp    ! because zu(nzb) is negative
2471                data_output_pr(i)     = data_output_pr(i)(2:)
2472             ENDIF
2473
2474          CASE ( 'w"vpt"' )
2475             dopr_index(i) = 45
2476             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2477             hom(:,2,45,:) = SPREAD( zw, 2, statistic_regions+1 )
2478
2479          CASE ( 'w*vpt*' )
2480             dopr_index(i) = 46
2481             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2482             hom(:,2,46,:) = SPREAD( zw, 2, statistic_regions+1 )
2483
2484          CASE ( 'wvpt' )
2485             dopr_index(i) = 47
2486             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2487             hom(:,2,47,:) = SPREAD( zw, 2, statistic_regions+1 )
2488
2489          CASE ( 'w"q"' )
2490             IF (  .NOT.  humidity )  THEN
2491                message_string = 'data_output_pr = ' //                        &
2492                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2493                                 'lemented for humidity = .FALSE.'
2494                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2495             ELSE
2496                dopr_index(i) = 48
2497                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2498                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2499             ENDIF
2500
2501          CASE ( 'w*q*' )
2502             IF (  .NOT.  humidity )  THEN
2503                message_string = 'data_output_pr = ' //                        &
2504                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2505                                 'lemented for humidity = .FALSE.'
2506                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2507             ELSE
2508                dopr_index(i) = 49
2509                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2510                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2511             ENDIF
2512
2513          CASE ( 'wq' )
2514             IF (  .NOT.  humidity )  THEN
2515                message_string = 'data_output_pr = ' //                        &
2516                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2517                                 'lemented for humidity = .FALSE.'
2518                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2519             ELSE
2520                dopr_index(i) = 50
2521                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2522                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2523             ENDIF
2524
2525          CASE ( 'w"s"' )
2526             IF (  .NOT.  passive_scalar )  THEN
2527                message_string = 'data_output_pr = ' //                        &
2528                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2529                                 'lemented for passive_scalar = .FALSE.'
2530                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2531             ELSE
2532                dopr_index(i) = 119
2533                dopr_unit(i)  = 'kg/m3 m/s'
2534                hom(:,2,119,:) = SPREAD( zw, 2, statistic_regions+1 )
2535             ENDIF
2536
2537          CASE ( 'w*s*' )
2538             IF (  .NOT.  passive_scalar )  THEN
2539                message_string = 'data_output_pr = ' //                        &
2540                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2541                                 'lemented for passive_scalar = .FALSE.'
2542                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2543             ELSE
2544                dopr_index(i) = 116
2545                dopr_unit(i)  = 'kg/m3 m/s'
2546                hom(:,2,116,:) = SPREAD( zw, 2, statistic_regions+1 )
2547             ENDIF
2548
2549          CASE ( 'ws' )
2550             IF (  .NOT.  passive_scalar )  THEN
2551                message_string = 'data_output_pr = ' //                        &
2552                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2553                                 'lemented for passive_scalar = .FALSE.'
2554                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2555             ELSE
2556                dopr_index(i) = 120
2557                dopr_unit(i)  = 'kg/m3 m/s'
2558                hom(:,2,120,:) = SPREAD( zw, 2, statistic_regions+1 )
2559             ENDIF
2560
2561          CASE ( 'w"qv"' )
2562             IF ( humidity  .AND.  .NOT.  cloud_physics )  THEN
2563                dopr_index(i) = 48
2564                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2565                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2566             ELSEIF ( humidity  .AND.  cloud_physics )  THEN
2567                dopr_index(i) = 51
2568                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2569                hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
2570             ELSE
2571                message_string = 'data_output_pr = ' //                        &
2572                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2573                                 'lemented for cloud_physics = .FALSE. an&' // &
2574                                 'd humidity = .FALSE.'
2575                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2576             ENDIF
2577
2578          CASE ( 'w*qv*' )
2579             IF ( humidity  .AND.  .NOT. cloud_physics )                       &
2580             THEN
2581                dopr_index(i) = 49
2582                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2583                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2584             ELSEIF( humidity .AND. cloud_physics ) THEN
2585                dopr_index(i) = 52
2586                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2587                hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
2588             ELSE
2589                message_string = 'data_output_pr = ' //                        &
2590                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2591                                 'lemented for cloud_physics = .FALSE. an&' // &
2592                                 'd humidity = .FALSE.'
2593                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2594             ENDIF
2595
2596          CASE ( 'wqv' )
2597             IF ( humidity  .AND.  .NOT.  cloud_physics )  THEN
2598                dopr_index(i) = 50
2599                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2600                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2601             ELSEIF ( humidity  .AND.  cloud_physics )  THEN
2602                dopr_index(i) = 53
2603                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2604                hom(:,2,53,:) = SPREAD( zw, 2, statistic_regions+1 )
2605             ELSE
2606                message_string = 'data_output_pr = ' //                        &
2607                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2608                                 'lemented for cloud_physics = .FALSE. an&' // &
2609                                 'd humidity = .FALSE.'
2610                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2611             ENDIF
2612
2613          CASE ( 'ql' )
2614             IF (  .NOT.  cloud_physics  .AND.  .NOT.  cloud_droplets )  THEN
2615                message_string = 'data_output_pr = ' //                        &
2616                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2617                                 'lemented for cloud_physics = .FALSE. or'  // &
2618                                 '&cloud_droplets = .FALSE.'
2619                CALL message( 'check_parameters', 'PA0096', 1, 2, 0, 6, 0 )
2620             ELSE
2621                dopr_index(i) = 54
2622                dopr_unit(i)  = 'kg/kg'
2623                hom(:,2,54,:)  = SPREAD( zu, 2, statistic_regions+1 )
2624             ENDIF
2625
2626          CASE ( 'w*u*u*:dz' )
2627             dopr_index(i) = 55
2628             dopr_unit(i)  = 'm2/s3'
2629             hom(:,2,55,:) = SPREAD( zu, 2, statistic_regions+1 )
2630
2631          CASE ( 'w*p*:dz' )
2632             dopr_index(i) = 56
2633             dopr_unit(i)  = 'm2/s3'
2634             hom(:,2,56,:) = SPREAD( zw, 2, statistic_regions+1 )
2635
2636          CASE ( 'w"e:dz' )
2637             dopr_index(i) = 57
2638             dopr_unit(i)  = 'm2/s3'
2639             hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 )
2640
2641
2642          CASE ( 'u"pt"' )
2643             dopr_index(i) = 58
2644             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2645             hom(:,2,58,:) = SPREAD( zu, 2, statistic_regions+1 )
2646
2647          CASE ( 'u*pt*' )
2648             dopr_index(i) = 59
2649             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2650             hom(:,2,59,:) = SPREAD( zu, 2, statistic_regions+1 )
2651
2652          CASE ( 'upt_t' )
2653             dopr_index(i) = 60
2654             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2655             hom(:,2,60,:) = SPREAD( zu, 2, statistic_regions+1 )
2656
2657          CASE ( 'v"pt"' )
2658             dopr_index(i) = 61
2659             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2660             hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
2661
2662          CASE ( 'v*pt*' )
2663             dopr_index(i) = 62
2664             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2665             hom(:,2,62,:) = SPREAD( zu, 2, statistic_regions+1 )
2666
2667          CASE ( 'vpt_t' )
2668             dopr_index(i) = 63
2669             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2670             hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
2671
2672          CASE ( 'rho_ocean' )
2673             IF (  .NOT.  ocean ) THEN
2674                message_string = 'data_output_pr = ' //                        &
2675                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2676                                 'lemented for ocean = .FALSE.'
2677                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2678             ELSE
2679                dopr_index(i) = 64
2680                dopr_unit(i)  = 'kg/m3'
2681                hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
2682                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2683                   dopr_initial_index(i) = 77
2684                   hom(:,2,77,:)         = SPREAD( zu, 2, statistic_regions+1 )
2685                   hom(nzb,2,77,:)       = 0.0_wp    ! because zu(nzb) is negative
2686                   data_output_pr(i)     = data_output_pr(i)(2:)
2687                ENDIF
2688             ENDIF
2689
2690          CASE ( 'w"sa"' )
2691             IF (  .NOT.  ocean ) THEN
2692                message_string = 'data_output_pr = ' //                        &
2693                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2694                                 'lemented for ocean = .FALSE.'
2695                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2696             ELSE
2697                dopr_index(i) = 65
2698                dopr_unit(i)  = 'psu m/s'
2699                hom(:,2,65,:) = SPREAD( zw, 2, statistic_regions+1 )
2700             ENDIF
2701
2702          CASE ( 'w*sa*' )
2703             IF (  .NOT. ocean  ) THEN
2704                message_string = 'data_output_pr = ' //                        &
2705                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2706                                 'lemented for ocean = .FALSE.'
2707                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2708             ELSE
2709                dopr_index(i) = 66
2710                dopr_unit(i)  = 'psu m/s'
2711                hom(:,2,66,:) = SPREAD( zw, 2, statistic_regions+1 )
2712             ENDIF
2713
2714          CASE ( 'wsa' )
2715             IF (  .NOT.  ocean ) THEN
2716                message_string = 'data_output_pr = ' //                        &
2717                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2718                                 'lemented for ocean = .FALSE.'
2719                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2720             ELSE
2721                dopr_index(i) = 67
2722                dopr_unit(i)  = 'psu m/s'
2723                hom(:,2,67,:) = SPREAD( zw, 2, statistic_regions+1 )
2724             ENDIF
2725
2726          CASE ( 'w*p*' )
2727             dopr_index(i) = 68
2728             dopr_unit(i)  = 'm3/s3'
2729             hom(:,2,68,:) = SPREAD( zu, 2, statistic_regions+1 )
2730
2731          CASE ( 'w"e' )
2732             dopr_index(i) = 69
2733             dopr_unit(i)  = 'm3/s3'
2734             hom(:,2,69,:) = SPREAD( zu, 2, statistic_regions+1 )
2735
2736          CASE ( 'q*2' )
2737             IF (  .NOT.  humidity )  THEN
2738                message_string = 'data_output_pr = ' //                        &
2739                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2740                                 'lemented for humidity = .FALSE.'
2741                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2742             ELSE
2743                dopr_index(i) = 70
2744                dopr_unit(i)  = 'kg2/kg2'
2745                hom(:,2,70,:) = SPREAD( zu, 2, statistic_regions+1 )
2746             ENDIF
2747
2748          CASE ( 'prho' )
2749             IF (  .NOT.  ocean ) THEN
2750                message_string = 'data_output_pr = ' //                        &
2751                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2752                                 'lemented for ocean = .FALSE.'
2753                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2754             ELSE
2755                dopr_index(i) = 71
2756                dopr_unit(i)  = 'kg/m3'
2757                hom(:,2,71,:) = SPREAD( zu, 2, statistic_regions+1 )
2758             ENDIF
2759
2760          CASE ( 'hyp' )
2761             dopr_index(i) = 72
2762             dopr_unit(i)  = 'hPa'
2763             hom(:,2,72,:) = SPREAD( zu, 2, statistic_regions+1 )
2764
2765          CASE ( 'rho_air' )
2766             dopr_index(i)  = 119
2767             dopr_unit(i)   = 'kg/m3'
2768             hom(:,2,119,:) = SPREAD( zu, 2, statistic_regions+1 )
2769
2770          CASE ( 'rho_air_zw' )
2771             dopr_index(i)  = 120
2772             dopr_unit(i)   = 'kg/m3'
2773             hom(:,2,120,:) = SPREAD( zw, 2, statistic_regions+1 )
2774
2775          CASE ( 'nc' )
2776             IF (  .NOT.  cloud_physics )  THEN
2777                message_string = 'data_output_pr = ' //                        &
2778                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2779                                 'lemented for cloud_physics = .FALSE.'
2780                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2781             ELSEIF ( .NOT.  microphysics_morrison )  THEN
2782                message_string = 'data_output_pr = ' //                        &
2783                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2784                                 'lemented for cloud_scheme /= morrison'
2785                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2786             ELSE
2787                dopr_index(i) = 89
2788                dopr_unit(i)  = '1/m3'
2789                hom(:,2,89,:)  = SPREAD( zu, 2, statistic_regions+1 )
2790             ENDIF
2791
2792          CASE ( 'nr' )
2793             IF (  .NOT.  cloud_physics )  THEN
2794                message_string = 'data_output_pr = ' //                        &
2795                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2796                                 'lemented for cloud_physics = .FALSE.'
2797                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2798             ELSEIF ( .NOT.  microphysics_seifert )  THEN
2799                message_string = 'data_output_pr = ' //                        &
2800                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2801                                 'lemented for cloud_scheme /= seifert_beheng'
2802                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2803             ELSE
2804                dopr_index(i) = 73
2805                dopr_unit(i)  = '1/m3'
2806                hom(:,2,73,:)  = SPREAD( zu, 2, statistic_regions+1 )
2807             ENDIF
2808
2809          CASE ( 'qr' )
2810             IF (  .NOT.  cloud_physics )  THEN
2811                message_string = 'data_output_pr = ' //                        &
2812                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2813                                 'lemented for cloud_physics = .FALSE.'
2814                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2815             ELSEIF ( .NOT.  microphysics_seifert )  THEN
2816                message_string = 'data_output_pr = ' //                        &
2817                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2818                                 'lemented for cloud_scheme /= seifert_beheng'
2819                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2820             ELSE
2821                dopr_index(i) = 74
2822                dopr_unit(i)  = 'kg/kg'
2823                hom(:,2,74,:)  = SPREAD( zu, 2, statistic_regions+1 )
2824             ENDIF
2825
2826          CASE ( 'qc' )
2827             IF (  .NOT.  cloud_physics )  THEN
2828                message_string = 'data_output_pr = ' //                        &
2829                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2830                                 'lemented for cloud_physics = .FALSE.'
2831                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2832             ELSE
2833                dopr_index(i) = 75
2834                dopr_unit(i)  = 'kg/kg'
2835                hom(:,2,75,:)  = SPREAD( zu, 2, statistic_regions+1 )
2836             ENDIF
2837
2838          CASE ( 'prr' )
2839             IF (  .NOT.  cloud_physics )  THEN
2840                message_string = 'data_output_pr = ' //                        &
2841                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2842                                 'lemented for cloud_physics = .FALSE.'
2843                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2844             ELSEIF ( microphysics_sat_adjust )  THEN
2845                message_string = 'data_output_pr = ' //                        &
2846                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
2847                                 'ilable for cloud_scheme = saturation_adjust'
2848                CALL message( 'check_parameters', 'PA0422', 1, 2, 0, 6, 0 )
2849             ELSE
2850                dopr_index(i) = 76
2851                dopr_unit(i)  = 'kg/kg m/s'
2852                hom(:,2,76,:)  = SPREAD( zu, 2, statistic_regions+1 )
2853             ENDIF
2854
2855          CASE ( 'ug' )
2856             dopr_index(i) = 78
2857             dopr_unit(i)  = 'm/s'
2858             hom(:,2,78,:) = SPREAD( zu, 2, statistic_regions+1 )
2859
2860          CASE ( 'vg' )
2861             dopr_index(i) = 79
2862             dopr_unit(i)  = 'm/s'
2863             hom(:,2,79,:) = SPREAD( zu, 2, statistic_regions+1 )
2864
2865          CASE ( 'w_subs' )
2866             IF (  .NOT.  large_scale_subsidence )  THEN
2867                message_string = 'data_output_pr = ' //                        &
2868                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2869                                 'lemented for large_scale_subsidence = .FALSE.'
2870                CALL message( 'check_parameters', 'PA0382', 1, 2, 0, 6, 0 )
2871             ELSE
2872                dopr_index(i) = 80
2873                dopr_unit(i)  = 'm/s'
2874                hom(:,2,80,:) = SPREAD( zu, 2, statistic_regions+1 )
2875             ENDIF
2876
2877          CASE ( 's*2' )
2878             IF (  .NOT.  passive_scalar )  THEN
2879                message_string = 'data_output_pr = ' //                        &
2880                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2881                                 'lemented for passive_scalar = .FALSE.'
2882                CALL message( 'check_parameters', 'PA0185', 1, 2, 0, 6, 0 )
2883             ELSE
2884                dopr_index(i) = 118
2885                dopr_unit(i)  = 'kg2/kg2'
2886                hom(:,2,118,:) = SPREAD( zu, 2, statistic_regions+1 )
2887             ENDIF
2888
2889
2890
2891          CASE DEFAULT
2892
2893             CALL lsm_check_data_output_pr( data_output_pr(i), i, unit,        &
2894                                            dopr_unit(i) )
2895
2896             CALL lsf_nudging_check_data_output_pr( data_output_pr(i), i,      &
2897                                                    unit, dopr_unit(i) )
2898
2899             IF ( unit == 'illegal' )  THEN
2900                CALL radiation_check_data_output_pr( data_output_pr(i), i,     &
2901                                                     unit, dopr_unit(i) )
2902             ENDIF
2903
2904             IF ( unit == 'illegal' )  THEN
2905                unit = ''
2906                CALL user_check_data_output_pr( data_output_pr(i), i, unit )
2907             ENDIF
2908
2909             IF ( unit == 'illegal' )  THEN
2910                IF ( data_output_pr_user(1) /= ' ' )  THEN
2911                   message_string = 'illegal value for data_output_pr or ' //  &
2912                                    'data_output_pr_user = "' //               &
2913                                    TRIM( data_output_pr(i) ) // '"'
2914                   CALL message( 'check_parameters', 'PA0097', 1, 2, 0, 6, 0 )
2915                ELSE
2916                   message_string = 'illegal value for data_output_pr = "' //  &
2917                                    TRIM( data_output_pr(i) ) // '"'
2918                   CALL message( 'check_parameters', 'PA0098', 1, 2, 0, 6, 0 )
2919                ENDIF
2920             ENDIF
2921
2922       END SELECT
2923
2924    ENDDO
2925
2926
2927!
2928!-- Append user-defined data output variables to the standard data output
2929    IF ( data_output_user(1) /= ' ' )  THEN
2930       i = 1
2931       DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 500 )
2932          i = i + 1
2933       ENDDO
2934       j = 1
2935       DO  WHILE ( data_output_user(j) /= ' '  .AND.  j <= 500 )
2936          IF ( i > 500 )  THEN
2937             message_string = 'number of output quantitities given by data' // &
2938                '_output and data_output_user exceeds the limit of 500'
2939             CALL message( 'check_parameters', 'PA0102', 1, 2, 0, 6, 0 )
2940          ENDIF
2941          data_output(i) = data_output_user(j)
2942          i = i + 1
2943          j = j + 1
2944       ENDDO
2945    ENDIF
2946
2947!
2948!-- Check and set steering parameters for 2d/3d data output and averaging
2949    i   = 1
2950    DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 500 )
2951!
2952!--    Check for data averaging
2953       ilen = LEN_TRIM( data_output(i) )
2954       j = 0                                                 ! no data averaging
2955       IF ( ilen > 3 )  THEN
2956          IF ( data_output(i)(ilen-2:ilen) == '_av' )  THEN
2957             j = 1                                           ! data averaging
2958             data_output(i) = data_output(i)(1:ilen-3)
2959          ENDIF
2960       ENDIF
2961!
2962!--    Check for cross section or volume data
2963       ilen = LEN_TRIM( data_output(i) )
2964       k = 0                                                   ! 3d data
2965       var = data_output(i)(1:ilen)
2966       IF ( ilen > 3 )  THEN
2967          IF ( data_output(i)(ilen-2:ilen) == '_xy'  .OR.                      &
2968               data_output(i)(ilen-2:ilen) == '_xz'  .OR.                      &
2969               data_output(i)(ilen-2:ilen) == '_yz' )  THEN
2970             k = 1                                             ! 2d data
2971             var = data_output(i)(1:ilen-3)
2972          ENDIF
2973       ENDIF
2974
2975!
2976!--    Check for allowed value and set units
2977       SELECT CASE ( TRIM( var ) )
2978
2979          CASE ( 'e' )
2980             IF ( constant_diffusion )  THEN
2981                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2982                                 'res constant_diffusion = .FALSE.'
2983                CALL message( 'check_parameters', 'PA0103', 1, 2, 0, 6, 0 )
2984             ENDIF
2985             unit = 'm2/s2'
2986
2987          CASE ( 'lpt' )
2988             IF (  .NOT.  cloud_physics )  THEN
2989                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2990                         'res cloud_physics = .TRUE.'
2991                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2992             ENDIF
2993             unit = 'K'
2994
2995          CASE ( 'nc' )
2996             IF (  .NOT.  cloud_physics )  THEN
2997                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2998                         'res cloud_physics = .TRUE.'
2999                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3000             ELSEIF ( .NOT.  microphysics_morrison )  THEN
3001                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3002                         'res = microphysics morrison '
3003                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
3004             ENDIF
3005             unit = '1/m3'
3006
3007          CASE ( 'nr' )
3008             IF (  .NOT.  cloud_physics )  THEN
3009                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3010                         'res cloud_physics = .TRUE.'
3011                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3012             ELSEIF ( .NOT.  microphysics_seifert )  THEN
3013                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3014                         'res cloud_scheme = seifert_beheng'
3015                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
3016             ENDIF
3017             unit = '1/m3'
3018
3019          CASE ( 'pc', 'pr' )
3020             IF (  .NOT.  particle_advection )  THEN
3021                message_string = 'output of "' // TRIM( var ) // '" requir' // &
3022                   'es a "particles_par"-NAMELIST in the parameter file (PARIN)'
3023                CALL message( 'check_parameters', 'PA0104', 1, 2, 0, 6, 0 )
3024             ENDIF
3025             IF ( TRIM( var ) == 'pc' )  unit = 'number'
3026             IF ( TRIM( var ) == 'pr' )  unit = 'm'
3027
3028          CASE ( 'prr' )
3029             IF (  .NOT.  cloud_physics )  THEN
3030                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3031                         'res cloud_physics = .TRUE.'
3032                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3033             ELSEIF ( microphysics_sat_adjust )  THEN
3034                message_string = 'output of "' // TRIM( var ) // '" is ' //  &
3035                         'not available for cloud_scheme = saturation_adjust'
3036                CALL message( 'check_parameters', 'PA0423', 1, 2, 0, 6, 0 )
3037             ENDIF
3038             unit = 'kg/kg m/s'
3039
3040          CASE ( 'q', 'vpt' )
3041             IF (  .NOT.  humidity )  THEN
3042                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3043                                 'res humidity = .TRUE.'
3044                CALL message( 'check_parameters', 'PA0105', 1, 2, 0, 6, 0 )
3045             ENDIF
3046             IF ( TRIM( var ) == 'q'   )  unit = 'kg/kg'
3047             IF ( TRIM( var ) == 'vpt' )  unit = 'K'
3048
3049          CASE ( 'qc' )
3050             IF (  .NOT.  cloud_physics )  THEN
3051                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3052                         'res cloud_physics = .TRUE.'
3053                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3054             ENDIF
3055             unit = 'kg/kg'
3056
3057          CASE ( 'ql' )
3058             IF ( .NOT.  ( cloud_physics  .OR.  cloud_droplets ) )  THEN
3059                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3060                         'res cloud_physics = .TRUE. or cloud_droplets = .TRUE.'
3061                CALL message( 'check_parameters', 'PA0106', 1, 2, 0, 6, 0 )
3062             ENDIF
3063             unit = 'kg/kg'
3064
3065          CASE ( 'ql_c', 'ql_v', 'ql_vp' )
3066             IF (  .NOT.  cloud_droplets )  THEN
3067                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3068                                 'res cloud_droplets = .TRUE.'
3069                CALL message( 'check_parameters', 'PA0107', 1, 2, 0, 6, 0 )
3070             ENDIF
3071             IF ( TRIM( var ) == 'ql_c'  )  unit = 'kg/kg'
3072             IF ( TRIM( var ) == 'ql_v'  )  unit = 'm3'
3073             IF ( TRIM( var ) == 'ql_vp' )  unit = 'none'
3074
3075          CASE ( 'qr' )
3076             IF (  .NOT.  cloud_physics )  THEN
3077                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3078                         'res cloud_physics = .TRUE.'
3079                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3080             ELSEIF ( .NOT.  microphysics_seifert ) THEN
3081                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3082                         'res cloud_scheme = seifert_beheng'
3083                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
3084             ENDIF
3085             unit = 'kg/kg'
3086
3087          CASE ( 'qv' )
3088             IF (  .NOT.  cloud_physics )  THEN
3089                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3090                                 'res cloud_physics = .TRUE.'
3091                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3092             ENDIF
3093             unit = 'kg/kg'
3094
3095          CASE ( 'rho_ocean' )
3096             IF (  .NOT.  ocean )  THEN
3097                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3098                                 'res ocean = .TRUE.'
3099                CALL message( 'check_parameters', 'PA0109', 1, 2, 0, 6, 0 )
3100             ENDIF
3101             unit = 'kg/m3'
3102
3103          CASE ( 's' )
3104             IF (  .NOT.  passive_scalar )  THEN
3105                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3106                                 'res passive_scalar = .TRUE.'
3107                CALL message( 'check_parameters', 'PA0110', 1, 2, 0, 6, 0 )
3108             ENDIF
3109             unit = 'conc'
3110
3111          CASE ( 'sa' )
3112             IF (  .NOT.  ocean )  THEN
3113                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3114                                 'res ocean = .TRUE.'
3115                CALL message( 'check_parameters', 'PA0109', 1, 2, 0, 6, 0 )
3116             ENDIF
3117             unit = 'psu'
3118
3119          CASE ( 'p', 'pt', 'u', 'v', 'w' )
3120             IF ( TRIM( var ) == 'p'  )  unit = 'Pa'
3121             IF ( TRIM( var ) == 'pt' )  unit = 'K'
3122             IF ( TRIM( var ) == 'u'  )  unit = 'm/s'
3123             IF ( TRIM( var ) == 'v'  )  unit = 'm/s'
3124             IF ( TRIM( var ) == 'w'  )  unit = 'm/s'
3125             CONTINUE
3126
3127          CASE ( 'lwp*', 'ol*', 'pra*', 'prr*', 'qsws*', 'shf*', 'ssws*', 't*', &
3128                 'u*', 'z0*', 'z0h*', 'z0q*' )
3129             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
3130                message_string = 'illegal value for data_output: "' //         &
3131                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
3132                                 'cross sections are allowed for this value'
3133                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
3134             ENDIF
3135
3136             IF ( TRIM( var ) == 'lwp*'  .AND.  .NOT. cloud_physics )  THEN
3137                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3138                                 'res cloud_physics = .TRUE.'
3139                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3140             ENDIF
3141             IF ( TRIM( var ) == 'pra*'  .AND.                                 &
3142                  .NOT. ( microphysics_kessler .OR. microphysics_seifert ) )  THEN
3143                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3144                                 'res cloud_scheme = kessler or seifert_beheng'
3145                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
3146             ENDIF
3147             IF ( TRIM( var ) == 'pra*'  .AND.  j == 1 )  THEN
3148                message_string = 'temporal averaging of precipitation ' //     &
3149                          'amount "' // TRIM( var ) // '" is not possible'
3150                CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 )
3151             ENDIF
3152             IF ( TRIM( var ) == 'prr*'  .AND.                                 &
3153                  .NOT. ( microphysics_kessler .OR. microphysics_seifert ) )  THEN
3154                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3155                                 'res cloud_scheme = kessler or seifert_beheng'
3156                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
3157             ENDIF
3158             IF ( TRIM( var ) == 'qsws*'  .AND.  .NOT.  humidity )  THEN
3159                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3160                                 'res humidity = .TRUE.'
3161                CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 )
3162             ENDIF
3163             IF ( TRIM( var ) == 'ssws*'  .AND.  .NOT.  passive_scalar )  THEN
3164                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3165                                 'res passive_scalar = .TRUE.'
3166                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
3167             ENDIF
3168
3169             IF ( TRIM( var ) == 'lwp*'   )  unit = 'kg/m2'
3170             IF ( TRIM( var ) == 'ol*'   )   unit = 'm'
3171             IF ( TRIM( var ) == 'pra*'   )  unit = 'mm'
3172             IF ( TRIM( var ) == 'prr*'   )  unit = 'mm/s'
3173             IF ( TRIM( var ) == 'qsws*'  )  unit = 'kgm/kgs'
3174             IF ( TRIM( var ) == 'shf*'   )  unit = 'K*m/s'
3175             IF ( TRIM( var ) == 'ssws*'  )  unit = 'kg/m2*s'
3176             IF ( TRIM( var ) == 't*'     )  unit = 'K'
3177             IF ( TRIM( var ) == 'u*'     )  unit = 'm/s'
3178             IF ( TRIM( var ) == 'z0*'    )  unit = 'm'
3179             IF ( TRIM( var ) == 'z0h*'   )  unit = 'm'
3180
3181          CASE DEFAULT
3182
3183             CALL lsm_check_data_output ( var, unit, i, ilen, k )
3184
3185             IF ( unit == 'illegal' )  THEN
3186                CALL radiation_check_data_output( var, unit, i, ilen, k )
3187             ENDIF
3188
3189!
3190!--          Block of urban surface model outputs
3191             IF ( unit == 'illegal' .AND. urban_surface .AND. var(1:4) == 'usm_' ) THEN
3192                 CALL usm_check_data_output( var, unit )
3193             ENDIF
3194
3195!
3196!--          Block of plant canopy model outputs
3197             IF ( unit == 'illegal' .AND. plant_canopy .AND. var(1:4) == 'pcm_' ) THEN
3198                 CALL pcm_check_data_output( var, unit )
3199             ENDIF
3200
3201             IF ( unit == 'illegal' )  THEN
3202                unit = ''
3203                CALL user_check_data_output( var, unit )
3204             ENDIF
3205
3206             IF ( unit == 'illegal' )  THEN
3207                IF ( data_output_user(1) /= ' ' )  THEN
3208                   message_string = 'illegal value for data_output or ' //     &
3209                         'data_output_user = "' // TRIM( data_output(i) ) // '"'
3210                   CALL message( 'check_parameters', 'PA0114', 1, 2, 0, 6, 0 )
3211                ELSE
3212                   message_string = 'illegal value for data_output = "' //     &
3213                                    TRIM( data_output(i) ) // '"'
3214                   CALL message( 'check_parameters', 'PA0115', 1, 2, 0, 6, 0 )
3215                ENDIF
3216             ENDIF
3217
3218       END SELECT
3219!
3220!--    Set the internal steering parameters appropriately
3221       IF ( k == 0 )  THEN
3222          do3d_no(j)              = do3d_no(j) + 1
3223          do3d(j,do3d_no(j))      = data_output(i)
3224          do3d_unit(j,do3d_no(j)) = unit
3225       ELSE
3226          do2d_no(j)              = do2d_no(j) + 1
3227          do2d(j,do2d_no(j))      = data_output(i)
3228          do2d_unit(j,do2d_no(j)) = unit
3229          IF ( data_output(i)(ilen-2:ilen) == '_xy' )  THEN
3230             data_output_xy(j) = .TRUE.
3231          ENDIF
3232          IF ( data_output(i)(ilen-2:ilen) == '_xz' )  THEN
3233             data_output_xz(j) = .TRUE.
3234          ENDIF
3235          IF ( data_output(i)(ilen-2:ilen) == '_yz' )  THEN
3236             data_output_yz(j) = .TRUE.
3237          ENDIF
3238       ENDIF
3239
3240       IF ( j == 1 )  THEN
3241!
3242!--       Check, if variable is already subject to averaging
3243          found = .FALSE.
3244          DO  k = 1, doav_n
3245             IF ( TRIM( doav(k) ) == TRIM( var ) )  found = .TRUE.
3246          ENDDO
3247
3248          IF ( .NOT. found )  THEN
3249             doav_n = doav_n + 1
3250             doav(doav_n) = var
3251          ENDIF
3252       ENDIF
3253
3254       i = i + 1
3255    ENDDO
3256
3257!
3258!-- Averaged 2d or 3d output requires that an averaging interval has been set
3259    IF ( doav_n > 0  .AND.  averaging_interval == 0.0_wp )  THEN
3260       WRITE( message_string, * )  'output of averaged quantity "',            &
3261                                   TRIM( doav(1) ), '_av" requires to set a ', &
3262                                   'non-zero & averaging interval'
3263       CALL message( 'check_parameters', 'PA0323', 1, 2, 0, 6, 0 )
3264    ENDIF
3265
3266!
3267!-- Check sectional planes and store them in one shared array
3268    IF ( ANY( section_xy > nz + 1 ) )  THEN
3269       WRITE( message_string, * )  'section_xy must be <= nz + 1 = ', nz + 1
3270       CALL message( 'check_parameters', 'PA0319', 1, 2, 0, 6, 0 )
3271    ENDIF
3272    IF ( ANY( section_xz > ny + 1 ) )  THEN
3273       WRITE( message_string, * )  'section_xz must be <= ny + 1 = ', ny + 1
3274       CALL message( 'check_parameters', 'PA0320', 1, 2, 0, 6, 0 )
3275    ENDIF
3276    IF ( ANY( section_yz > nx + 1 ) )  THEN
3277       WRITE( message_string, * )  'section_yz must be <= nx + 1 = ', nx + 1
3278       CALL message( 'check_parameters', 'PA0321', 1, 2, 0, 6, 0 )
3279    ENDIF
3280    section(:,1) = section_xy
3281    section(:,2) = section_xz
3282    section(:,3) = section_yz
3283
3284!
3285!-- Upper plot limit for 2D vertical sections
3286    IF ( z_max_do2d == -1.0_wp )  z_max_do2d = zu(nzt)
3287    IF ( z_max_do2d < zu(nzb+1)  .OR.  z_max_do2d > zu(nzt) )  THEN
3288       WRITE( message_string, * )  'z_max_do2d = ', z_max_do2d,                &
3289                    ' must be >= ', zu(nzb+1), '(zu(nzb+1)) and <= ', zu(nzt), &
3290                    ' (zu(nzt))'
3291       CALL message( 'check_parameters', 'PA0116', 1, 2, 0, 6, 0 )
3292    ENDIF
3293
3294!
3295!-- Upper plot limit for 3D arrays
3296    IF ( nz_do3d == -9999 )  nz_do3d = nzt + 1
3297
3298!
3299!-- Set output format string (used in header)
3300    SELECT CASE ( netcdf_data_format )
3301       CASE ( 1 )
3302          netcdf_data_format_string = 'netCDF classic'
3303       CASE ( 2 )
3304          netcdf_data_format_string = 'netCDF 64bit offset'
3305       CASE ( 3 )
3306          netcdf_data_format_string = 'netCDF4/HDF5'
3307       CASE ( 4 )
3308          netcdf_data_format_string = 'netCDF4/HDF5 classic'
3309       CASE ( 5 )
3310          netcdf_data_format_string = 'parallel netCDF4/HDF5'
3311       CASE ( 6 )
3312          netcdf_data_format_string = 'parallel netCDF4/HDF5 classic'
3313
3314    END SELECT
3315
3316!
3317!-- Check mask conditions
3318    DO mid = 1, max_masks
3319       IF ( data_output_masks(mid,1) /= ' '  .OR.                              &
3320            data_output_masks_user(mid,1) /= ' ' ) THEN
3321          masks = masks + 1
3322       ENDIF
3323    ENDDO
3324
3325    IF ( masks < 0  .OR.  masks > max_masks )  THEN
3326       WRITE( message_string, * )  'illegal value: masks must be >= 0 and ',   &
3327            '<= ', max_masks, ' (=max_masks)'
3328       CALL message( 'check_parameters', 'PA0325', 1, 2, 0, 6, 0 )
3329    ENDIF
3330    IF ( masks > 0 )  THEN
3331       mask_scale(1) = mask_scale_x
3332       mask_scale(2) = mask_scale_y
3333       mask_scale(3) = mask_scale_z
3334       IF ( ANY( mask_scale <= 0.0_wp ) )  THEN
3335          WRITE( message_string, * )                                           &
3336               'illegal value: mask_scale_x, mask_scale_y and mask_scale_z',   &
3337               'must be > 0.0'
3338          CALL message( 'check_parameters', 'PA0326', 1, 2, 0, 6, 0 )
3339       ENDIF
3340!
3341!--    Generate masks for masked data output
3342!--    Parallel netcdf output is not tested so far for masked data, hence
3343!--    netcdf_data_format is switched back to non-paralell output.
3344       netcdf_data_format_save = netcdf_data_format
3345       IF ( netcdf_data_format > 4 )  THEN
3346          IF ( netcdf_data_format == 5 ) netcdf_data_format = 3
3347          IF ( netcdf_data_format == 6 ) netcdf_data_format = 4
3348          message_string = 'netCDF file formats '//                            &
3349                           '5 (parallel netCDF 4) and ' //                     &
3350                           '6 (parallel netCDF 4 Classic model) '//            &
3351                           '&are currently not supported (not yet tested) ' // &
3352                           'for masked data.&Using respective non-parallel' // &
3353                           ' output for masked data.'
3354          CALL message( 'check_parameters', 'PA0383', 0, 0, 0, 6, 0 )
3355       ENDIF
3356       CALL init_masks
3357       netcdf_data_format = netcdf_data_format_save
3358    ENDIF
3359
3360!
3361!-- Check the NetCDF data format
3362    IF ( netcdf_data_format > 2 )  THEN
3363#if defined( __netcdf4 )
3364       CONTINUE
3365#else
3366       message_string = 'netCDF: netCDF4 format requested but no ' //          &
3367                        'cpp-directive __netcdf4 given & switch '  //          &
3368                        'back to 64-bit offset format'
3369       CALL message( 'check_parameters', 'PA0171', 0, 1, 0, 6, 0 )
3370       netcdf_data_format = 2
3371#endif
3372    ENDIF
3373    IF ( netcdf_data_format > 4 )  THEN
3374#if defined( __netcdf4 ) && defined( __netcdf4_parallel )
3375       CONTINUE
3376#else
3377       message_string = 'netCDF: netCDF4 parallel output requested but no ' // &
3378                        'cpp-directive __netcdf4_parallel given & switch '  // &
3379                        'back to netCDF4 non-parallel output'
3380       CALL message( 'check_parameters', 'PA0099', 0, 1, 0, 6, 0 )
3381       netcdf_data_format = netcdf_data_format - 2
3382#endif
3383    ENDIF
3384
3385!
3386!-- Calculate fixed number of output time levels for parallel netcdf output.
3387!-- The time dimension has to be defined as limited for parallel output,
3388!-- because otherwise the I/O performance drops significantly.
3389    IF ( netcdf_data_format > 4 )  THEN
3390
3391!
3392!--    Check if any of the follwoing data output interval is 0.0s, which is
3393!--    not allowed for parallel output.
3394       CALL check_dt_do( dt_do3d,    'dt_do3d'    )
3395       CALL check_dt_do( dt_do2d_xy, 'dt_do2d_xy' )
3396       CALL check_dt_do( dt_do2d_xz, 'dt_do2d_xz' )
3397       CALL check_dt_do( dt_do2d_yz, 'dt_do2d_yz' )
3398
3399       ntdim_3d(0)    = INT( ( end_time - skip_time_do3d ) / dt_do3d )
3400       IF ( do3d_at_begin ) ntdim_3d(0) = ntdim_3d(0) + 1
3401       ntdim_3d(1)    = INT( ( end_time - skip_time_data_output_av )           &
3402                             / dt_data_output_av )
3403       ntdim_2d_xy(0) = INT( ( end_time - skip_time_do2d_xy ) / dt_do2d_xy )
3404       ntdim_2d_xz(0) = INT( ( end_time - skip_time_do2d_xz ) / dt_do2d_xz )
3405       ntdim_2d_yz(0) = INT( ( end_time - skip_time_do2d_yz ) / dt_do2d_yz )
3406       IF ( do2d_at_begin )  THEN
3407          ntdim_2d_xy(0) = ntdim_2d_xy(0) + 1
3408          ntdim_2d_xz(0) = ntdim_2d_xz(0) + 1
3409          ntdim_2d_yz(0) = ntdim_2d_yz(0) + 1
3410       ENDIF
3411       ntdim_2d_xy(1) = ntdim_3d(1)
3412       ntdim_2d_xz(1) = ntdim_3d(1)
3413       ntdim_2d_yz(1) = ntdim_3d(1)
3414
3415    ENDIF
3416
3417!
3418!-- Check, whether a constant diffusion coefficient shall be used
3419    IF ( km_constant /= -1.0_wp )  THEN
3420       IF ( km_constant < 0.0_wp )  THEN
3421          WRITE( message_string, * )  'km_constant = ', km_constant, ' < 0.0'
3422          CALL message( 'check_parameters', 'PA0121', 1, 2, 0, 6, 0 )
3423       ELSE
3424          IF ( prandtl_number < 0.0_wp )  THEN
3425             WRITE( message_string, * )  'prandtl_number = ', prandtl_number,  &
3426                                         ' < 0.0'
3427             CALL message( 'check_parameters', 'PA0122', 1, 2, 0, 6, 0 )
3428          ENDIF
3429          constant_diffusion = .TRUE.
3430
3431          IF ( constant_flux_layer )  THEN
3432             message_string = 'constant_flux_layer is not allowed with fixed ' &
3433                              // 'value of km'
3434             CALL message( 'check_parameters', 'PA0123', 1, 2, 0, 6, 0 )
3435          ENDIF
3436       ENDIF
3437    ENDIF
3438
3439!
3440!-- In case of non-cyclic lateral boundaries and a damping layer for the
3441!-- potential temperature, check the width of the damping layer
3442    IF ( bc_lr /= 'cyclic' ) THEN
3443       IF ( pt_damping_width < 0.0_wp  .OR.                                    &
3444            pt_damping_width > REAL( (nx+1) * dx ) )  THEN
3445          message_string = 'pt_damping_width out of range'
3446          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
3447       ENDIF
3448    ENDIF
3449
3450    IF ( bc_ns /= 'cyclic' )  THEN
3451       IF ( pt_damping_width < 0.0_wp  .OR.                                    &
3452            pt_damping_width > REAL( (ny+1) * dy ) )  THEN
3453          message_string = 'pt_damping_width out of range'
3454          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
3455       ENDIF
3456    ENDIF
3457
3458!
3459!-- Check value range for zeta = z/L
3460    IF ( zeta_min >= zeta_max )  THEN
3461       WRITE( message_string, * )  'zeta_min = ', zeta_min, ' must be less ',  &
3462                                   'than zeta_max = ', zeta_max
3463       CALL message( 'check_parameters', 'PA0125', 1, 2, 0, 6, 0 )
3464    ENDIF
3465
3466!
3467!-- Check random generator
3468    IF ( (random_generator /= 'system-specific'      .AND.                     &
3469          random_generator /= 'random-parallel'   )  .AND.                     &
3470          random_generator /= 'numerical-recipes' )  THEN
3471       message_string = 'unknown random generator: random_generator = "' //    &
3472                        TRIM( random_generator ) // '"'
3473       CALL message( 'check_parameters', 'PA0135', 1, 2, 0, 6, 0 )
3474    ENDIF
3475
3476!
3477!-- Determine upper and lower hight level indices for random perturbations
3478    IF ( disturbance_level_b == -9999999.9_wp )  THEN
3479       IF ( ocean )  THEN
3480          disturbance_level_b     = zu((nzt*2)/3)
3481          disturbance_level_ind_b = ( nzt * 2 ) / 3
3482       ELSE
3483          disturbance_level_b     = zu(nzb+3)
3484          disturbance_level_ind_b = nzb + 3
3485       ENDIF
3486    ELSEIF ( disturbance_level_b < zu(3) )  THEN
3487       WRITE( message_string, * )  'disturbance_level_b = ',                   &
3488                           disturbance_level_b, ' must be >= ', zu(3), '(zu(3))'
3489       CALL message( 'check_parameters', 'PA0126', 1, 2, 0, 6, 0 )
3490    ELSEIF ( disturbance_level_b > zu(nzt-2) )  THEN
3491       WRITE( message_string, * )  'disturbance_level_b = ',                   &
3492                   disturbance_level_b, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
3493       CALL message( 'check_parameters', 'PA0127', 1, 2, 0, 6, 0 )
3494    ELSE
3495       DO  k = 3, nzt-2
3496          IF ( disturbance_level_b <= zu(k) )  THEN
3497             disturbance_level_ind_b = k
3498             EXIT
3499          ENDIF
3500       ENDDO
3501    ENDIF
3502
3503    IF ( disturbance_level_t == -9999999.9_wp )  THEN
3504       IF ( ocean )  THEN
3505          disturbance_level_t     = zu(nzt-3)
3506          disturbance_level_ind_t = nzt - 3
3507       ELSE
3508          disturbance_level_t     = zu(nzt/3)
3509          disturbance_level_ind_t = nzt / 3
3510       ENDIF
3511    ELSEIF ( disturbance_level_t > zu(nzt-2) )  THEN
3512       WRITE( message_string, * )  'disturbance_level_t = ',                   &
3513                   disturbance_level_t, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
3514       CALL message( 'check_parameters', 'PA0128', 1, 2, 0, 6, 0 )
3515    ELSEIF ( disturbance_level_t < disturbance_level_b )  THEN
3516       WRITE( message_string, * )  'disturbance_level_t = ',                   &
3517                   disturbance_level_t, ' must be >= disturbance_level_b = ',  &
3518                   disturbance_level_b
3519       CALL message( 'check_parameters', 'PA0129', 1, 2, 0, 6, 0 )
3520    ELSE
3521       DO  k = 3, nzt-2
3522          IF ( disturbance_level_t <= zu(k) )  THEN
3523             disturbance_level_ind_t = k
3524             EXIT
3525          ENDIF
3526       ENDDO
3527    ENDIF
3528
3529!
3530!-- Check again whether the levels determined this way are ok.
3531!-- Error may occur at automatic determination and too few grid points in
3532!-- z-direction.
3533    IF ( disturbance_level_ind_t < disturbance_level_ind_b )  THEN
3534       WRITE( message_string, * )  'disturbance_level_ind_t = ',               &
3535                disturbance_level_ind_t, ' must be >= disturbance_level_ind_b = ', &
3536                disturbance_level_ind_b
3537       CALL message( 'check_parameters', 'PA0130', 1, 2, 0, 6, 0 )
3538    ENDIF
3539
3540!
3541!-- Determine the horizontal index range for random perturbations.
3542!-- In case of non-cyclic horizontal boundaries, no perturbations are imposed
3543!-- near the inflow and the perturbation area is further limited to ...(1)
3544!-- after the initial phase of the flow.
3545
3546    IF ( bc_lr /= 'cyclic' )  THEN
3547       IF ( inflow_disturbance_begin == -1 )  THEN
3548          inflow_disturbance_begin = MIN( 10, nx/2 )
3549       ENDIF
3550       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > nx )&
3551       THEN
3552          message_string = 'inflow_disturbance_begin out of range'
3553          CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
3554       ENDIF
3555       IF ( inflow_disturbance_end == -1 )  THEN
3556          inflow_disturbance_end = MIN( 100, 3*nx/4 )
3557       ENDIF
3558       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > nx )    &
3559       THEN
3560          message_string = 'inflow_disturbance_end out of range'
3561          CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
3562       ENDIF
3563    ELSEIF ( bc_ns /= 'cyclic' )  THEN
3564       IF ( inflow_disturbance_begin == -1 )  THEN
3565          inflow_disturbance_begin = MIN( 10, ny/2 )
3566       ENDIF
3567       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > ny )&
3568       THEN
3569          message_string = 'inflow_disturbance_begin out of range'
3570          CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
3571       ENDIF
3572       IF ( inflow_disturbance_end == -1 )  THEN
3573          inflow_disturbance_end = MIN( 100, 3*ny/4 )
3574       ENDIF
3575       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > ny )    &
3576       THEN
3577          message_string = 'inflow_disturbance_end out of range'
3578          CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
3579       ENDIF
3580    ENDIF
3581
3582    IF ( random_generator == 'random-parallel' )  THEN
3583       dist_nxl = nxl;  dist_nxr = nxr
3584       dist_nys = nys;  dist_nyn = nyn
3585       IF ( bc_lr == 'radiation/dirichlet' )  THEN
3586          dist_nxr    = MIN( nx - inflow_disturbance_begin, nxr )
3587          dist_nxl(1) = MAX( nx - inflow_disturbance_end, nxl )
3588       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
3589          dist_nxl    = MAX( inflow_disturbance_begin, nxl )
3590          dist_nxr(1) = MIN( inflow_disturbance_end, nxr )
3591       ELSEIF ( bc_lr == 'nested' )  THEN
3592          dist_nxl = MAX( inflow_disturbance_begin, nxl )
3593       ENDIF
3594       IF ( bc_ns == 'dirichlet/radiation' )  THEN
3595          dist_nyn    = MIN( ny - inflow_disturbance_begin, nyn )
3596          dist_nys(1) = MAX( ny - inflow_disturbance_end, nys )
3597       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
3598          dist_nys    = MAX( inflow_disturbance_begin, nys )
3599          dist_nyn(1) = MIN( inflow_disturbance_end, nyn )
3600       ELSEIF ( bc_ns == 'nested' )  THEN
3601          dist_nys = MAX( inflow_disturbance_begin, nys )
3602       ENDIF
3603    ELSE
3604       dist_nxl = 0;  dist_nxr = nx
3605       dist_nys = 0;  dist_nyn = ny
3606       IF ( bc_lr == 'radiation/dirichlet' )  THEN
3607          dist_nxr    = nx - inflow_disturbance_begin
3608          dist_nxl(1) = nx - inflow_disturbance_end
3609       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
3610          dist_nxl    = inflow_disturbance_begin
3611          dist_nxr(1) = inflow_disturbance_end
3612       ENDIF
3613       IF ( bc_ns == 'dirichlet/radiation' )  THEN
3614          dist_nyn    = ny - inflow_disturbance_begin
3615          dist_nys(1) = ny - inflow_disturbance_end
3616       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
3617          dist_nys    = inflow_disturbance_begin
3618          dist_nyn(1) = inflow_disturbance_end
3619       ENDIF
3620    ENDIF
3621
3622!
3623!-- A turbulent inflow requires Dirichlet conditions at the respective inflow
3624!-- boundary (so far, a turbulent inflow is realized from the left side only)
3625    IF ( turbulent_inflow  .AND.  bc_lr /= 'dirichlet/radiation' )  THEN
3626       message_string = 'turbulent_inflow = .T. requires a Dirichlet ' //      &
3627                        'condition at the inflow boundary'
3628       CALL message( 'check_parameters', 'PA0133', 1, 2, 0, 6, 0 )
3629    ENDIF
3630
3631!
3632!-- Turbulent inflow requires that 3d arrays have been cyclically filled with
3633!-- data from prerun in the first main run
3634    IF ( turbulent_inflow  .AND.  initializing_actions /= 'cyclic_fill'  .AND. &
3635         initializing_actions /= 'read_restart_data' )  THEN
3636       message_string = 'turbulent_inflow = .T. requires ' //                  &
3637                        'initializing_actions = ''cyclic_fill'' or ' //        &
3638                        'initializing_actions = ''read_restart_data'' '
3639       CALL message( 'check_parameters', 'PA0055', 1, 2, 0, 6, 0 )
3640    ENDIF
3641
3642!
3643!-- In case of turbulent inflow calculate the index of the recycling plane
3644    IF ( turbulent_inflow )  THEN
3645       IF ( recycling_width <= dx  .OR.  recycling_width >= nx * dx )  THEN
3646          WRITE( message_string, * )  'illegal value for recycling_width:', &
3647                                      ' ', recycling_width
3648          CALL message( 'check_parameters', 'PA0134', 1, 2, 0, 6, 0 )
3649       ENDIF
3650!
3651!--    Calculate the index
3652       recycling_plane = recycling_width / dx
3653!
3654!--    Because the y-shift is done with a distance of INT( npey / 2 ) no shift
3655!--    is possible if there is only one PE in y direction.
3656       IF ( recycling_yshift .AND. pdims(2) < 2 )  THEN
3657          WRITE( message_string, * )  'recycling_yshift = .T. requires more',  &
3658                                      ' than one processor in y direction'
3659          CALL message( 'check_parameters', 'PA0421', 1, 2, 0, 6, 0 )
3660       ENDIF
3661    ENDIF
3662
3663
3664    IF ( turbulent_outflow )  THEN
3665!
3666!--    Turbulent outflow requires Dirichlet conditions at the respective inflow
3667!--    boundary (so far, a turbulent outflow is realized at the right side only)
3668       IF ( bc_lr /= 'dirichlet/radiation' )  THEN
3669          message_string = 'turbulent_outflow = .T. requires ' //              &
3670                           'bc_lr = "dirichlet/radiation"'
3671          CALL message( 'check_parameters', 'PA0038', 1, 2, 0, 6, 0 )
3672       ENDIF
3673!
3674!--    The ouflow-source plane must lay inside the model domain
3675       IF ( outflow_source_plane < dx  .OR.  &
3676            outflow_source_plane > nx * dx )  THEN
3677          WRITE( message_string, * )  'illegal value for outflow_source'//     &
3678                                      '_plane: ', outflow_source_plane
3679          CALL message( 'check_parameters', 'PA0145', 1, 2, 0, 6, 0 )
3680       ENDIF
3681    ENDIF
3682
3683!
3684!-- Determine damping level index for 1D model
3685    IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
3686       IF ( damp_level_1d == -1.0_wp )  THEN
3687          damp_level_1d     = zu(nzt+1)
3688          damp_level_ind_1d = nzt + 1
3689       ELSEIF ( damp_level_1d < 0.0_wp  .OR.  damp_level_1d > zu(nzt+1) )  THEN
3690          WRITE( message_string, * )  'damp_level_1d = ', damp_level_1d,       &
3691                 ' must be >= 0.0 and <= ', zu(nzt+1), '(zu(nzt+1))'
3692          CALL message( 'check_parameters', 'PA0136', 1, 2, 0, 6, 0 )
3693       ELSE
3694          DO  k = 1, nzt+1
3695             IF ( damp_level_1d <= zu(k) )  THEN
3696                damp_level_ind_1d = k
3697                EXIT
3698             ENDIF
3699          ENDDO
3700       ENDIF
3701    ENDIF
3702
3703!
3704!-- Check some other 1d-model parameters
3705    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND.                   &
3706         TRIM( mixing_length_1d ) /= 'blackadar' )  THEN
3707       message_string = 'mixing_length_1d = "' // TRIM( mixing_length_1d ) //  &
3708                        '" is unknown'
3709       CALL message( 'check_parameters', 'PA0137', 1, 2, 0, 6, 0 )
3710    ENDIF
3711    IF ( TRIM( dissipation_1d ) /= 'as_in_3d_model'  .AND.                     &
3712         TRIM( dissipation_1d ) /= 'detering' )  THEN
3713       message_string = 'dissipation_1d = "' // TRIM( dissipation_1d ) //      &
3714                        '" is unknown'
3715       CALL message( 'check_parameters', 'PA0138', 1, 2, 0, 6, 0 )
3716    ENDIF
3717
3718!
3719!-- Set time for the next user defined restart (time_restart is the
3720!-- internal parameter for steering restart events)
3721    IF ( restart_time /= 9999999.9_wp )  THEN
3722       IF ( restart_time > time_since_reference_point )  THEN
3723          time_restart = restart_time
3724       ENDIF
3725    ELSE
3726!
3727!--    In case of a restart run, set internal parameter to default (no restart)
3728!--    if the NAMELIST-parameter restart_time is omitted
3729       time_restart = 9999999.9_wp
3730    ENDIF
3731
3732!
3733!-- Check pressure gradient conditions
3734    IF ( dp_external  .AND.  conserve_volume_flow )  THEN
3735       WRITE( message_string, * )  'Both dp_external and conserve_volume_flo', &
3736            'w are .TRUE. but one of them must be .FALSE.'
3737       CALL message( 'check_parameters', 'PA0150', 1, 2, 0, 6, 0 )
3738    ENDIF
3739    IF ( dp_external )  THEN
3740       IF ( dp_level_b < zu(nzb)  .OR.  dp_level_b > zu(nzt) )  THEN
3741          WRITE( message_string, * )  'dp_level_b = ', dp_level_b, ' is out ', &
3742               ' of range'
3743          CALL message( 'check_parameters', 'PA0151', 1, 2, 0, 6, 0 )
3744       ENDIF
3745       IF ( .NOT. ANY( dpdxy /= 0.0_wp ) )  THEN
3746          WRITE( message_string, * )  'dp_external is .TRUE. but dpdxy is ze', &
3747               'ro, i.e. the external pressure gradient & will not be applied'
3748          CALL message( 'check_parameters', 'PA0152', 0, 1, 0, 6, 0 )
3749       ENDIF
3750    ENDIF
3751    IF ( ANY( dpdxy /= 0.0_wp )  .AND.  .NOT.  dp_external )  THEN
3752       WRITE( message_string, * )  'dpdxy is nonzero but dp_external is ',     &
3753            '.FALSE., i.e. the external pressure gradient & will not be applied'
3754       CALL message( 'check_parameters', 'PA0153', 0, 1, 0, 6, 0 )
3755    ENDIF
3756    IF ( conserve_volume_flow )  THEN
3757       IF ( TRIM( conserve_volume_flow_mode ) == 'default' )  THEN
3758
3759          conserve_volume_flow_mode = 'initial_profiles'
3760
3761       ELSEIF ( TRIM( conserve_volume_flow_mode ) /= 'initial_profiles' .AND.  &
3762            TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' )  THEN
3763          WRITE( message_string, * )  'unknown conserve_volume_flow_mode: ',   &
3764               conserve_volume_flow_mode
3765          CALL message( 'check_parameters', 'PA0154', 1, 2, 0, 6, 0 )
3766       ENDIF
3767       IF ( (bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic')  .AND.                &
3768          TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
3769          WRITE( message_string, * )  'non-cyclic boundary conditions ',       &
3770               'require  conserve_volume_flow_mode = ''initial_profiles'''
3771          CALL message( 'check_parameters', 'PA0155', 1, 2, 0, 6, 0 )
3772       ENDIF
3773    ENDIF
3774    IF ( ( u_bulk /= 0.0_wp  .OR.  v_bulk /= 0.0_wp )  .AND.                   &
3775         ( .NOT. conserve_volume_flow  .OR.                                    &
3776         TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' ) )  THEN
3777       WRITE( message_string, * )  'nonzero bulk velocity requires ',          &
3778            'conserve_volume_flow = .T. and ',                                 &
3779            'conserve_volume_flow_mode = ''bulk_velocity'''
3780       CALL message( 'check_parameters', 'PA0157', 1, 2, 0, 6, 0 )
3781    ENDIF
3782
3783!
3784!-- Check particle attributes
3785    IF ( particle_color /= 'none' )  THEN
3786       IF ( particle_color /= 'absuv'  .AND.  particle_color /= 'pt*'  .AND.   &
3787            particle_color /= 'z' )  THEN
3788          message_string = 'illegal value for parameter particle_color: ' //   &
3789                           TRIM( particle_color)
3790          CALL message( 'check_parameters', 'PA0313', 1, 2, 0, 6, 0 )
3791       ELSE
3792          IF ( color_interval(2) <= color_interval(1) )  THEN
3793             message_string = 'color_interval(2) <= color_interval(1)'
3794             CALL message( 'check_parameters', 'PA0315', 1, 2, 0, 6, 0 )
3795          ENDIF
3796       ENDIF
3797    ENDIF
3798
3799    IF ( particle_dvrpsize /= 'none' )  THEN
3800       IF ( particle_dvrpsize /= 'absw' )  THEN
3801          message_string = 'illegal value for parameter particle_dvrpsize:' // &
3802                           ' ' // TRIM( particle_color)
3803          CALL message( 'check_parameters', 'PA0314', 1, 2, 0, 6, 0 )
3804       ELSE
3805          IF ( dvrpsize_interval(2) <= dvrpsize_interval(1) )  THEN
3806             message_string = 'dvrpsize_interval(2) <= dvrpsize_interval(1)'
3807             CALL message( 'check_parameters', 'PA0316', 1, 2, 0, 6, 0 )
3808          ENDIF
3809       ENDIF
3810    ENDIF
3811
3812!
3813!-- Prevent empty time records in volume, cross-section and masked data in case
3814!-- of non-parallel netcdf-output in restart runs
3815    IF ( netcdf_data_format < 5 )  THEN
3816       IF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
3817          do3d_time_count    = 0
3818          do2d_xy_time_count = 0
3819          do2d_xz_time_count = 0
3820          do2d_yz_time_count = 0
3821          domask_time_count  = 0
3822       ENDIF
3823    ENDIF
3824
3825!
3826!-- Check for valid setting of most_method
3827    IF ( TRIM( most_method ) /= 'circular'  .AND.                              &
3828         TRIM( most_method ) /= 'newton'    .AND.                              &
3829         TRIM( most_method ) /= 'lookup' )  THEN
3830       message_string = 'most_method = "' // TRIM( most_method ) //            &
3831                        '" is unknown'
3832       CALL message( 'check_parameters', 'PA0416', 1, 2, 0, 6, 0 )
3833    ENDIF
3834
3835!
3836!-- Check roughness length, which has to be smaller than dz/2
3837    IF ( ( constant_flux_layer .OR.  &
3838           INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) &
3839         .AND. roughness_length >= 0.5 * dz )  THEN
3840       message_string = 'roughness_length must be smaller than dz/2'
3841       CALL message( 'check_parameters', 'PA0424', 1, 2, 0, 6, 0 )
3842    ENDIF
3843
3844    CALL location_message( 'finished', .TRUE. )
3845!
3846!-- Check &userpar parameters
3847    CALL user_check_parameters
3848
3849 CONTAINS
3850
3851!------------------------------------------------------------------------------!
3852! Description:
3853! ------------
3854!> Check the length of data output intervals. In case of parallel NetCDF output
3855!> the time levels of the output files need to be fixed. Therefore setting the
3856!> output interval to 0.0s (usually used to output each timestep) is not
3857!> possible as long as a non-fixed timestep is used.
3858!------------------------------------------------------------------------------!
3859
3860    SUBROUTINE check_dt_do( dt_do, dt_do_name )
3861
3862       IMPLICIT NONE
3863
3864       CHARACTER (LEN=*), INTENT (IN) :: dt_do_name !< parin variable name
3865
3866       REAL(wp), INTENT (INOUT)       :: dt_do      !< data output interval
3867
3868       IF ( dt_do == 0.0_wp )  THEN
3869          IF ( dt_fixed )  THEN
3870             WRITE( message_string, '(A,F9.4,A)' )  'Output at every '  //     &
3871                    'timestep is wanted (' // dt_do_name // ' = 0.0).&'//      &
3872                    'Setting the output interval to the fixed timestep '//     &
3873                    'dt = ', dt, 's.'
3874             CALL message( 'check_parameters', 'PA0060', 0, 0, 0, 6, 0 )
3875             dt_do = dt
3876          ELSE
3877             message_string = dt_do_name // ' = 0.0 while using a ' //         &
3878                              'variable timestep and parallel netCDF4 ' //     &
3879                              'is not allowed.'
3880             CALL message( 'check_parameters', 'PA0081', 1, 2, 0, 6, 0 )
3881          ENDIF
3882       ENDIF
3883
3884    END SUBROUTINE check_dt_do
3885
3886
3887
3888
3889!------------------------------------------------------------------------------!
3890! Description:
3891! ------------
3892!> Inititalizes the vertical profiles of scalar quantities.
3893!------------------------------------------------------------------------------!
3894
3895    SUBROUTINE init_vertical_profiles( vertical_gradient_level_ind,            &
3896                                       vertical_gradient_level,                &
3897                                       vertical_gradient,                      &
3898                                       pr_init, surface_value, bc_t_val )
3899
3900
3901       IMPLICIT NONE
3902
3903       INTEGER(iwp) ::  i     !< counter
3904       INTEGER(iwp), DIMENSION(1:10) ::  vertical_gradient_level_ind !< vertical grid indices for gradient levels
3905
3906       REAL(wp)     ::  bc_t_val      !< model top gradient
3907       REAL(wp)     ::  gradient      !< vertica gradient of the respective quantity
3908       REAL(wp)     ::  surface_value !< surface value of the respecitve quantity
3909
3910
3911       REAL(wp), DIMENSION(0:nz+1) ::  pr_init                 !< initialisation profile
3912       REAL(wp), DIMENSION(1:10)   ::  vertical_gradient       !< given vertical gradient
3913       REAL(wp), DIMENSION(1:10)   ::  vertical_gradient_level !< given vertical gradient level
3914
3915       i = 1
3916       gradient = 0.0_wp
3917
3918       IF (  .NOT.  ocean )  THEN
3919
3920          vertical_gradient_level_ind(1) = 0
3921          DO  k = 1, nzt+1
3922             IF ( i < 11 )  THEN
3923                IF ( vertical_gradient_level(i) < zu(k)  .AND.            &
3924                     vertical_gradient_level(i) >= 0.0_wp )  THEN
3925                   gradient = vertical_gradient(i) / 100.0_wp
3926                   vertical_gradient_level_ind(i) = k - 1
3927                   i = i + 1
3928                ENDIF
3929             ENDIF
3930             IF ( gradient /= 0.0_wp )  THEN
3931                IF ( k /= 1 )  THEN
3932                   pr_init(k) = pr_init(k-1) + dzu(k) * gradient
3933                ELSE
3934                   pr_init(k) = pr_init(k-1) + dzu(k) * gradient
3935                ENDIF
3936             ELSE
3937                pr_init(k) = pr_init(k-1)
3938             ENDIF
3939   !
3940   !--       Avoid negative values
3941             IF ( pr_init(k) < 0.0_wp )  THEN
3942                pr_init(k) = 0.0_wp
3943             ENDIF
3944          ENDDO
3945
3946       ELSE
3947
3948          vertical_gradient_level_ind(1) = nzt+1
3949          DO  k = nzt, 0, -1
3950             IF ( i < 11 )  THEN
3951                IF ( vertical_gradient_level(i) > zu(k)  .AND.            &
3952                     vertical_gradient_level(i) <= 0.0_wp )  THEN
3953                   gradient = vertical_gradient(i) / 100.0_wp
3954                   vertical_gradient_level_ind(i) = k + 1
3955                   i = i + 1
3956                ENDIF
3957             ENDIF
3958             IF ( gradient /= 0.0_wp )  THEN
3959                IF ( k /= nzt )  THEN
3960                   pr_init(k) = pr_init(k+1) - dzu(k+1) * gradient
3961                ELSE
3962                   pr_init(k)   = surface_value - 0.5_wp * dzu(k+1) * gradient
3963                   pr_init(k+1) = surface_value + 0.5_wp * dzu(k+1) * gradient
3964                ENDIF
3965             ELSE
3966                pr_init(k) = pr_init(k+1)
3967             ENDIF
3968   !
3969   !--       Avoid negative humidities
3970             IF ( pr_init(k) < 0.0_wp )  THEN
3971                pr_init(k) = 0.0_wp
3972             ENDIF
3973          ENDDO
3974
3975       ENDIF
3976
3977!
3978!--    In case of no given gradients, choose zero gradient conditions
3979       IF ( vertical_gradient_level(1) == -1.0_wp )  THEN
3980          vertical_gradient_level(1) = 0.0_wp
3981       ENDIF
3982!
3983!--    Store gradient at the top boundary for possible Neumann boundary condition
3984       bc_t_val  = ( pr_init(nzt+1) - pr_init(nzt) ) / dzu(nzt+1)
3985
3986    END SUBROUTINE init_vertical_profiles
3987
3988
3989
3990!------------------------------------------------------------------------------!
3991! Description:
3992! ------------
3993!> Set the bottom and top boundary conditions for humidity and scalars.
3994!------------------------------------------------------------------------------!
3995
3996    SUBROUTINE set_bc_scalars( sq, bc_b, bc_t, ibc_b, ibc_t, err_nr_b, err_nr_t )
3997
3998
3999       IMPLICIT NONE
4000
4001       CHARACTER (LEN=1)   ::  sq                       !<
4002       CHARACTER (LEN=*)   ::  bc_b
4003       CHARACTER (LEN=*)   ::  bc_t
4004       CHARACTER (LEN=*)   ::  err_nr_b
4005       CHARACTER (LEN=*)   ::  err_nr_t
4006
4007       INTEGER(iwp)        ::  ibc_b
4008       INTEGER(iwp)        ::  ibc_t
4009
4010!
4011!--    Set Integer flags and check for possilbe errorneous settings for bottom
4012!--    boundary condition
4013       IF ( bc_b == 'dirichlet' )  THEN
4014          ibc_b = 0
4015       ELSEIF ( bc_b == 'neumann' )  THEN
4016          ibc_b = 1
4017       ELSE
4018          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
4019                           '_b ="' // TRIM( bc_b ) // '"'
4020          CALL message( 'check_parameters', err_nr_b, 1, 2, 0, 6, 0 )
4021       ENDIF
4022!
4023!--    Set Integer flags and check for possilbe errorneous settings for top
4024!--    boundary condition
4025       IF ( bc_t == 'dirichlet' )  THEN
4026          ibc_t = 0
4027       ELSEIF ( bc_t == 'neumann' )  THEN
4028          ibc_t = 1
4029       ELSEIF ( bc_t == 'initial_gradient' )  THEN
4030          ibc_t = 2
4031       ELSEIF ( bc_t == 'nested' )  THEN
4032          ibc_t = 3
4033       ELSE
4034          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
4035                           '_t ="' // TRIM( bc_t ) // '"'
4036          CALL message( 'check_parameters', err_nr_t, 1, 2, 0, 6, 0 )
4037       ENDIF
4038
4039
4040    END SUBROUTINE set_bc_scalars
4041
4042
4043
4044!------------------------------------------------------------------------------!
4045! Description:
4046! ------------
4047!> Check for consistent settings of bottom boundary conditions for humidity
4048!> and scalars.
4049!------------------------------------------------------------------------------!
4050
4051    SUBROUTINE check_bc_scalars( sq, bc_b, ibc_b,                 &
4052                                 err_nr_1, err_nr_2,       &
4053                                 constant_flux, surface_initial_change )
4054
4055
4056       IMPLICIT NONE
4057
4058       CHARACTER (LEN=1)   ::  sq                       !<
4059       CHARACTER (LEN=*)   ::  bc_b
4060       CHARACTER (LEN=*)   ::  err_nr_1
4061       CHARACTER (LEN=*)   ::  err_nr_2
4062
4063       INTEGER(iwp)        ::  ibc_b
4064
4065       LOGICAL             ::  constant_flux
4066
4067       REAL(wp)            ::  surface_initial_change
4068
4069!
4070!--    A given surface value implies Dirichlet boundary condition for
4071!--    the respective quantity. In this case specification of a constant flux is
4072!--    forbidden. However, an exception is made for large-scale forcing as well
4073!--    as land-surface model.
4074       IF ( .NOT. land_surface  .AND.  .NOT. large_scale_forcing )  THEN
4075          IF ( ibc_b == 0  .AND.  constant_flux )  THEN
4076             message_string = 'boundary condition: bc_' // TRIM( sq ) // '_b ' //  &
4077                              '= "' // TRIM( bc_b ) // '" is not allowed with ' // &
4078                              'prescribed surface flux'
4079             CALL message( 'check_parameters', err_nr_1, 1, 2, 0, 6, 0 )
4080          ENDIF
4081       ENDIF
4082       IF ( constant_flux  .AND.  surface_initial_change /= 0.0_wp )  THEN
4083          WRITE( message_string, * )  'a prescribed surface flux is not allo', &
4084                 'wed with ', sq, '_surface_initial_change (/=0) = ',          &
4085                 surface_initial_change
4086          CALL message( 'check_parameters', err_nr_2, 1, 2, 0, 6, 0 )
4087       ENDIF
4088
4089
4090    END SUBROUTINE check_bc_scalars
4091
4092
4093
4094 END SUBROUTINE check_parameters
Note: See TracBrowser for help on using the repository browser.