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

Last change on this file since 3766 was 3766, checked in by raasch, 3 years ago

unused_variables removed, bugfix in im_define_netcdf_grid argument list, trim added to avoid truncation compiler warnings, save attribute added to local targets to avoid outlive pointer target warning, first argument removed from module_interface_rrd_*, file module_interface reformatted with respect to coding standards, bugfix in surface_data_output_rrd_local (variable k removed)

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