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

Last change on this file since 2039 was 2039, checked in by gronemeier, 8 years ago

Increased the number of possible statistic regions to 99

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