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

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

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