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

Last change on this file since 3748 was 3735, checked in by dom_dwd_user, 5 years ago

biometeorology_mod.f90:
(N) Fixed auto-setting of thermal index calculation flags by output as
originally proposed by resler.
(C) removed bio_pet and outher configuration variables.
(C) Updated namelist.
(B) Forcing initialization of tmrt_av_grid to avoid mysterious mrt
values at i==0, j==0

module_interface_mod.f90:
(C) Receiving parameter j (averaging 0==.F./1==.T.) in
module_interface_check_data_output from check_parameters.f90.
(C) Passing j to bio_check_parameters.

check_parameters.f90:
(C) Passing j to module_interface_check_data_output

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