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

Last change on this file since 2246 was 2232, checked in by suehring, 7 years ago

Adjustments according new topography and surface-modelling concept implemented

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