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

Last change on this file since 1918 was 1918, checked in by raasch, 5 years ago

bugfixes for calculating run control quantities, bugfix for calculating pressure with fft-method in case of Neumann conditions both at bottom and top, steering of pres modified, ocean mode now uses initial density profile as reference in the buoyancy term

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