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

Last change on this file since 4039 was 4039, checked in by suehring, 2 years ago

diagnostic output: Modularize diagnostic output, rename subroutines; formatting adjustments; allocate arrays only when required; add output of uu, vv, ww to enable variance calculation via temporal EC method; radiation: bugfix in masked data output; flow_statistics: Correct conversion to kinematic vertical scalar fluxes in case of pw-scheme and statistic regions

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