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

Last change on this file since 2041 was 2041, checked in by gronemeier, 5 years ago

added missing svn-id string

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