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

Last change on this file since 4048 was 4048, checked in by knoop, 2 years ago

Moved turbulence_closure_mod calls into module_interface

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