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

Last change on this file since 1817 was 1817, checked in by maronga, 9 years ago

further modularization of land surface model

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