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

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

Nesting for chemical species implemented; Bugfix passive scalar boundary condition after anterpolation; Timeseries output of surface temperature; Enable initialization of 3D topography (was commented out so far)

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