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

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

Output of ground-heat flux at natural- and urban-type surfaces in one output variable; enable restart data of _av variables that belong to both land- and urban-surface model

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