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

Last change on this file since 3638 was 3637, checked in by knoop, 5 years ago

M Makefile

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