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

Last change on this file since 1825 was 1825, checked in by gronemeier, 8 years ago

last commit documented

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