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

Last change on this file since 1762 was 1762, checked in by hellstea, 6 years ago

Introduction of nested domain system

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