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

Last change on this file since 1822 was 1822, checked in by hoffmann, 8 years ago

changes in LPM and bulk cloud microphysics

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