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

Last change on this file since 2742 was 2742, checked in by suehring, 7 years ago

Output of surface temperature tsurf* at urban- and natural-type surfaces

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