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

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

pmc-change in server-client get-put, spectra-directives removed, spectra-package modularized

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