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

Last change on this file since 1804 was 1804, checked in by maronga, 5 years ago

removed parameter file check. update of mrungui for compilation with qt5

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