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

Last change on this file since 2851 was 2851, checked in by maronga, 6 years ago

bugfix: calculation of output time levels for restart runs (parallel NetCDF)

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