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

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

last commit documented

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