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

Last change on this file since 2044 was 2044, checked in by knoop, 7 years ago

Added error code for anelastic approximation

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