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

Last change on this file since 2037 was 2037, checked in by knoop, 8 years ago

Anelastic approximation implemented

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