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

Last change on this file since 2932 was 2932, checked in by maronga, 4 years ago

renamed all Fortran NAMELISTS

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