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

Last change on this file since 2174 was 2170, checked in by suehring, 8 years ago

last commit documented

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