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

Last change on this file since 1701 was 1701, checked in by maronga, 8 years ago

minor bugfixes for radiation model. bugfix in subjob

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