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

Last change on this file since 3705 was 3705, checked in by suehring, 4 years ago

last commit documented

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