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

Last change on this file since 3885 was 3885, checked in by kanani, 2 years ago

restructure/add location/debug messages

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