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

Last change on this file since 3291 was 3274, checked in by knoop, 6 years ago

Modularization of all bulk cloud physics code components

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