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

Last change on this file since 1914 was 1914, checked in by witha, 8 years ago

Merged branch/forwind into trunk

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