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

Last change on this file since 2365 was 2365, checked in by kanani, 7 years ago

Vertical nesting implemented (SadiqHuq?)

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