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

Last change on this file since 1692 was 1692, checked in by maronga, 9 years ago

last commit documented

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