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

Last change on this file since 4173 was 4172, checked in by oliver.maas, 5 years ago

Added optional method for recycling of absolute values of pt and q instead of recycling only the deviations from the mean profile at the recycling plane. With the new method two problems are solved: 1. A horizontally homogeneous temperature and humidity field is achieved, even for non-neutral boundary layers and thus no thermal circulation is triggered. 2. No gravity waves build up at the inflow due to cyclic boundary conditions for pt and q.

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