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

Last change on this file since 2225 was 2210, checked in by kanani, 7 years ago

last commit documented

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