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

Last change on this file since 1834 was 1834, checked in by raasch, 8 years ago

last commit documented

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