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

Last change on this file since 1682 was 1682, checked in by knoop, 9 years ago

Code annotations made doxygen readable

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