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

Last change on this file since 1952 was 1932, checked in by suehring, 8 years ago

last commit documented

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