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

Last change on this file since 2723 was 2718, checked in by maronga, 7 years ago

deleting of deprecated files; headers updated where needed

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