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

Last change on this file since 2338 was 2338, checked in by gronemeier, 7 years ago

modularized 1d model

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