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

Last change on this file since 3467 was 3467, checked in by suehring, 3 years ago

Branch salsa @3446 re-integrated into trunk

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