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

Last change on this file since 2186 was 2179, checked in by hellstea, 8 years ago

last commit documented

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