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

Last change on this file since 2737 was 2735, checked in by suehring, 6 years ago

Output of resistance also urban-type surfaces

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