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

Last change on this file since 2169 was 2169, checked in by suehring, 5 years ago

Bugfix, move setting of topography_grid_convention to init_grid, else, in case of restarts generic topography may not set properly if no value is prescribed

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