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

Last change on this file since 3464 was 3454, checked in by schwenkel, 5 years ago

Bugfix: Missing air_chemistry statement for profile check

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