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

Last change on this file since 1789 was 1789, checked in by maronga, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 184.5 KB
Line 
1!> @file check_parameters.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2015 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: check_parameters.f90 1789 2016-03-10 11:02:40Z maronga $
26!
27! 1788 2016-03-10 11:01:04Z maronga
28! Added check for use of most_method = 'lookup' in combination with water
29! surface presribed in the land surface model. Added output of z0q.
30! Syntax layout improved.
31!
32! 1786 2016-03-08 05:49:27Z raasch
33! cpp-direktives for spectra removed, check of spectra level removed
34!
35! 1783 2016-03-06 18:36:17Z raasch
36! netcdf variables and module name changed,
37! check of netcdf precision removed (is done in the netcdf module)
38!
39! 1764 2016-02-28 12:45:19Z raasch
40! output of nest id in run description header,
41! bugfix: check of validity of lateral boundary conditions moved to parin
42!
43! 1762 2016-02-25 12:31:13Z hellstea
44! Introduction of nested domain feature
45!
46! 1745 2016-02-05 13:06:51Z gronemeier
47! Bugfix: check data output intervals to be /= 0.0 in case of parallel NetCDF4
48!
49! 1701 2015-11-02 07:43:04Z maronga
50! Bugfix: definition of rad_net timeseries was missing
51!
52! 1691 2015-10-26 16:17:44Z maronga
53! Added output of Obukhov length (ol) and radiative heating rates for RRTMG.
54! Added checks for use of radiation / lsm with topography.
55!
56! 1682 2015-10-07 23:56:08Z knoop
57! Code annotations made doxygen readable
58!
59! 1606 2015-06-29 10:43:37Z maronga
60! Added check for use of RRTMG without netCDF.
61!
62! 1587 2015-05-04 14:19:01Z maronga
63! Added check for set of albedos when using RRTMG
64!
65! 1585 2015-04-30 07:05:52Z maronga
66! Added support for RRTMG
67!
68! 1575 2015-03-27 09:56:27Z raasch
69! multigrid_fast added as allowed pressure solver
70!
71! 1573 2015-03-23 17:53:37Z suehring
72! Bugfix: check for advection schemes in case of non-cyclic boundary conditions
73!
74! 1557 2015-03-05 16:43:04Z suehring
75! Added checks for monotonic limiter
76!
77! 1555 2015-03-04 17:44:27Z maronga
78! Added output of r_a and r_s. Renumbering of LSM PA-messages.
79!
80! 1553 2015-03-03 17:33:54Z maronga
81! Removed check for missing soil temperature profile as default values were added.
82!
83! 1551 2015-03-03 14:18:16Z maronga
84! Added various checks for land surface and radiation model. In the course of this
85! action, the length of the variable var has to be increased
86!
87! 1504 2014-12-04 09:23:49Z maronga
88! Bugfix: check for large_scale forcing got lost.
89!
90! 1500 2014-12-03 17:42:41Z maronga
91! Boundary conditions changed to dirichlet for land surface model
92!
93! 1496 2014-12-02 17:25:50Z maronga
94! Added checks for the land surface model
95!
96! 1484 2014-10-21 10:53:05Z kanani
97! Changes due to new module structure of the plant canopy model:
98!   module plant_canopy_model_mod added,
99!   checks regarding usage of new method for leaf area density profile
100!   construction added,
101!   lad-profile construction moved to new subroutine init_plant_canopy within
102!   the module plant_canopy_model_mod,
103!   drag_coefficient renamed to canopy_drag_coeff.
104! Missing KIND-attribute for REAL constant added
105!
106! 1455 2014-08-29 10:47:47Z heinze
107! empty time records in volume, cross-section and masked data output prevented 
108! in case of non-parallel netcdf-output in restart runs
109!
110! 1429 2014-07-15 12:53:45Z knoop
111! run_description_header exended to provide ensemble_member_nr if specified
112!
113! 1425 2014-07-05 10:57:53Z knoop
114! bugfix: perturbation domain modified for parallel random number generator
115!
116! 1402 2014-05-09 14:25:13Z raasch
117! location messages modified
118!
119! 1400 2014-05-09 14:03:54Z knoop
120! Check random generator extended by option random-parallel
121!
122! 1384 2014-05-02 14:31:06Z raasch
123! location messages added
124!
125! 1365 2014-04-22 15:03:56Z boeske
126! Usage of large scale forcing for pt and q enabled:
127! output for profiles of large scale advection (td_lsa_lpt, td_lsa_q),
128! large scale subsidence (td_sub_lpt, td_sub_q)
129! and nudging tendencies (td_nud_lpt, td_nud_q, td_nud_u and td_nud_v) added,
130! check if use_subsidence_tendencies is used correctly
131!
132! 1361 2014-04-16 15:17:48Z hoffmann
133! PA0363 removed
134! PA0362 changed
135!
136! 1359 2014-04-11 17:15:14Z hoffmann
137! Do not allow the execution of PALM with use_particle_tails, since particle
138! tails are currently not supported by our new particle structure.
139!
140! PA0084 not necessary for new particle structure
141!
142! 1353 2014-04-08 15:21:23Z heinze
143! REAL constants provided with KIND-attribute
144!
145! 1330 2014-03-24 17:29:32Z suehring
146! In case of SGS-particle velocity advection of TKE is also allowed with
147! dissipative 5th-order scheme.
148!
149! 1327 2014-03-21 11:00:16Z raasch
150! "baroclinicity" renamed "baroclinity", "ocean version" replaced by "ocean mode"
151! bugfix: duplicate error message 56 removed,
152! check of data_output_format and do3d_compress removed
153!
154! 1322 2014-03-20 16:38:49Z raasch
155! some REAL constants defined as wp-kind
156!
157! 1320 2014-03-20 08:40:49Z raasch
158! Kind-parameters added to all INTEGER and REAL declaration statements,
159! kinds are defined in new module kinds,
160! revision history before 2012 removed,
161! comment fields (!:) to be used for variable explanations added to
162! all variable declaration statements
163!
164! 1308 2014-03-13 14:58:42Z fricke
165! +netcdf_data_format_save
166! Calculate fixed number of output time levels for parallel netcdf output.
167! For masked data, parallel netcdf output is not tested so far, hence
168! netcdf_data_format is switched back to non-paralell output.
169!
170! 1299 2014-03-06 13:15:21Z heinze
171! enable usage of large_scale subsidence in combination with large_scale_forcing
172! output for profile of large scale vertical velocity w_subs added
173!
174! 1276 2014-01-15 13:40:41Z heinze
175! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
176!
177! 1241 2013-10-30 11:36:58Z heinze
178! output for profiles of ug and vg added
179! set surface_heatflux and surface_waterflux also in dependence on
180! large_scale_forcing
181! checks for nudging and large scale forcing from external file
182!
183! 1236 2013-09-27 07:21:13Z raasch
184! check number of spectra levels
185!
186! 1216 2013-08-26 09:31:42Z raasch
187! check for transpose_compute_overlap (temporary)
188!
189! 1214 2013-08-21 12:29:17Z kanani
190! additional check for simultaneous use of vertical grid stretching
191! and particle advection
192!
193! 1212 2013-08-15 08:46:27Z raasch
194! checks for poisfft_hybrid removed
195!
196! 1210 2013-08-14 10:58:20Z raasch
197! check for fftw
198!
199! 1179 2013-06-14 05:57:58Z raasch
200! checks and settings of buoyancy parameters and switches revised,
201! initial profile for rho added to hom (id=77)
202!
203! 1174 2013-05-31 10:28:08Z gryschka
204! Bugfix in computing initial profiles for ug, vg, lad, q in case of Atmosphere
205!
206! 1159 2013-05-21 11:58:22Z fricke
207! bc_lr/ns_dirneu/neudir removed
208!
209! 1115 2013-03-26 18:16:16Z hoffmann
210! unused variables removed
211! drizzle can be used without precipitation
212!
213! 1111 2013-03-08 23:54:10Z raasch
214! ibc_p_b = 2 removed
215!
216! 1103 2013-02-20 02:15:53Z raasch
217! Bugfix: turbulent inflow must not require cyclic fill in restart runs
218!
219! 1092 2013-02-02 11:24:22Z raasch
220! unused variables removed
221!
222! 1069 2012-11-28 16:18:43Z maronga
223! allow usage of topography in combination with cloud physics
224!
225! 1065 2012-11-22 17:42:36Z hoffmann
226! Bugfix: It is not allowed to use cloud_scheme = seifert_beheng without
227!         precipitation in order to save computational resources.
228!
229! 1060 2012-11-21 07:19:51Z raasch
230! additional check for parameter turbulent_inflow
231!
232! 1053 2012-11-13 17:11:03Z hoffmann
233! necessary changes for the new two-moment cloud physics scheme added:
234! - check cloud physics scheme (Kessler or Seifert and Beheng)
235! - plant_canopy is not allowed
236! - currently, only cache loop_optimization is allowed
237! - initial profiles of nr, qr
238! - boundary condition of nr, qr
239! - check output quantities (qr, nr, prr)
240!
241! 1036 2012-10-22 13:43:42Z raasch
242! code put under GPL (PALM 3.9)
243!
244! 1031/1034 2012-10-22 11:32:49Z raasch
245! check of netcdf4 parallel file support
246!
247! 1019 2012-09-28 06:46:45Z raasch
248! non-optimized version of prognostic_equations not allowed any more
249!
250! 1015 2012-09-27 09:23:24Z raasch
251! acc allowed for loop optimization,
252! checks for adjustment of mixing length to the Prandtl mixing length removed
253!
254! 1003 2012-09-14 14:35:53Z raasch
255! checks for cases with unequal subdomain sizes removed
256!
257! 1001 2012-09-13 14:08:46Z raasch
258! all actions concerning leapfrog- and upstream-spline-scheme removed
259!
260! 996 2012-09-07 10:41:47Z raasch
261! little reformatting
262
263! 978 2012-08-09 08:28:32Z fricke
264! setting of bc_lr/ns_dirneu/neudir
265! outflow damping layer removed
266! check for z0h*
267! check for pt_damping_width
268!
269! 964 2012-07-26 09:14:24Z raasch
270! check of old profil-parameters removed
271!
272! 940 2012-07-09 14:31:00Z raasch
273! checks for parameter neutral
274!
275! 924 2012-06-06 07:44:41Z maronga
276! Bugfix: preprocessor directives caused error during compilation
277!
278! 892 2012-05-02 13:51:44Z maronga
279! Bugfix for parameter file check ( excluding __netcdf4 )
280!
281! 866 2012-03-28 06:44:41Z raasch
282! use only 60% of the geostrophic wind as translation speed in case of Galilean
283! transformation and use_ug_for_galilei_tr = .T. in order to mimimize the
284! timestep
285!
286! 861 2012-03-26 14:18:34Z suehring
287! Check for topography and ws-scheme removed.
288! Check for loop_optimization = 'vector' and ws-scheme removed.
289!
290! 845 2012-03-07 10:23:05Z maronga
291! Bugfix: exclude __netcdf4 directive part from namelist file check compilation
292!
293! 828 2012-02-21 12:00:36Z raasch
294! check of collision_kernel extended
295!
296! 825 2012-02-19 03:03:44Z raasch
297! check for collision_kernel and curvature_solution_effects
298!
299! 809 2012-01-30 13:32:58Z maronga
300! Bugfix: replaced .AND. and .NOT. with && and ! in the preprocessor directives
301!
302! 807 2012-01-25 11:53:51Z maronga
303! New cpp directive "__check" implemented which is used by check_namelist_files
304!
305! Revision 1.1  1997/08/26 06:29:23  raasch
306! Initial revision
307!
308!
309! Description:
310! ------------
311!> Check control parameters and deduce further quantities.
312!------------------------------------------------------------------------------!
313 SUBROUTINE check_parameters
314 
315
316    USE arrays_3d
317    USE cloud_parameters
318    USE constants
319    USE control_parameters
320    USE dvrp_variables
321    USE grid_variables
322    USE indices
323    USE land_surface_model_mod
324    USE kinds
325    USE model_1d
326    USE netcdf_interface,                                                      &
327        ONLY:  dopr_unit, do2d_unit, do3d_unit, netcdf_data_format,            &
328               netcdf_data_format_string
329    USE particle_attributes
330    USE pegrid
331    USE plant_canopy_model_mod
332    USE pmc_interface,                                                         &
333        ONLY:  cpl_id, nested_run
334    USE profil_parameter
335    USE radiation_model_mod
336    USE statistics
337    USE subsidence_mod
338    USE statistics
339    USE transpose_indices
340
341
342    IMPLICIT NONE
343
344    CHARACTER (LEN=1)   ::  sq                       !<
345    CHARACTER (LEN=15)  ::  var                      !<
346    CHARACTER (LEN=7)   ::  unit                     !<
347    CHARACTER (LEN=8)   ::  date                     !<
348    CHARACTER (LEN=10)  ::  time                     !<
349    CHARACTER (LEN=10)  ::  ensemble_string          !<
350    CHARACTER (LEN=15)  ::  nest_string              !<
351    CHARACTER (LEN=40)  ::  coupling_string          !<
352    CHARACTER (LEN=100) ::  action                   !<
353
354    INTEGER(iwp) ::  i                               !<
355    INTEGER(iwp) ::  ilen                            !<
356    INTEGER(iwp) ::  iremote = 0                     !<
357    INTEGER(iwp) ::  j                               !<
358    INTEGER(iwp) ::  k                               !<
359    INTEGER(iwp) ::  kk                              !<
360    INTEGER(iwp) ::  netcdf_data_format_save         !<
361    INTEGER(iwp) ::  position                        !<
362    INTEGER(iwp) ::  prec                            !<
363   
364    LOGICAL     ::  found                            !<
365   
366    REAL(wp)    ::  gradient                         !<
367    REAL(wp)    ::  remote = 0.0_wp                  !<
368    REAL(wp)    ::  simulation_time_since_reference  !<
369
370
371    CALL location_message( 'checking parameters', .FALSE. )
372
373!
374!-- Check for overlap combinations, which are not realized yet
375    IF ( transpose_compute_overlap )  THEN
376       IF ( numprocs == 1 )  STOP '+++ transpose-compute-overlap not implemented for single PE runs'
377#if defined( __openacc )
378       STOP '+++ transpose-compute-overlap not implemented for GPU usage'
379#endif
380    ENDIF
381
382!
383!-- Warning, if host is not set
384    IF ( host(1:1) == ' ' )  THEN
385       message_string = '"host" is not set. Please check that environment ' // &
386                        'variable "localhost" & is set before running PALM'
387       CALL message( 'check_parameters', 'PA0001', 0, 0, 0, 6, 0 )
388    ENDIF
389
390!
391!-- Check the coupling mode
392    IF ( coupling_mode /= 'uncoupled'            .AND.  &
393         coupling_mode /= 'atmosphere_to_ocean'  .AND.  &
394         coupling_mode /= 'ocean_to_atmosphere' )  THEN
395       message_string = 'illegal coupling mode: ' // TRIM( coupling_mode )
396       CALL message( 'check_parameters', 'PA0002', 1, 2, 0, 6, 0 )
397    ENDIF
398
399!
400!-- Check dt_coupling, restart_time, dt_restart, end_time, dx, dy, nx and ny
401    IF ( coupling_mode /= 'uncoupled')  THEN
402
403       IF ( dt_coupling == 9999999.9_wp )  THEN
404          message_string = 'dt_coupling is not set but required for coup' //   &
405                           'ling mode "' //  TRIM( coupling_mode ) // '"'
406          CALL message( 'check_parameters', 'PA0003', 1, 2, 0, 6, 0 )
407       ENDIF
408
409#if defined( __parallel )
410
411!
412!--    NOTE: coupled runs have not been implemented in the check_namelist_files
413!--    program.
414!--    check_namelist_files will need the following information of the other
415!--    model (atmosphere/ocean).
416!       dt_coupling = remote
417!       dt_max = remote
418!       restart_time = remote
419!       dt_restart= remote
420!       simulation_time_since_reference = remote
421!       dx = remote
422
423
424#if ! defined( __check )
425       IF ( myid == 0 ) THEN
426          CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter,  &
427                         ierr )
428          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter,       &
429                         status, ierr )
430       ENDIF
431       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
432#endif     
433       IF ( dt_coupling /= remote )  THEN
434          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
435                 '": dt_coupling = ', dt_coupling, '& is not equal to ',       &
436                 'dt_coupling_remote = ', remote
437          CALL message( 'check_parameters', 'PA0004', 1, 2, 0, 6, 0 )
438       ENDIF
439       IF ( dt_coupling <= 0.0_wp )  THEN
440#if ! defined( __check )
441          IF ( myid == 0  ) THEN
442             CALL MPI_SEND( dt_max, 1, MPI_REAL, target_id, 19, comm_inter, ierr )
443             CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter,    &
444                            status, ierr )
445          ENDIF   
446          CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
447#endif         
448          dt_coupling = MAX( dt_max, remote )
449          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
450                 '": dt_coupling <= 0.0 & is not allowed and is reset to ',    &
451                 'MAX(dt_max(A,O)) = ', dt_coupling
452          CALL message( 'check_parameters', 'PA0005', 0, 1, 0, 6, 0 )
453       ENDIF
454#if ! defined( __check )
455       IF ( myid == 0 ) THEN
456          CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, &
457                         ierr )
458          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter,       &
459                         status, ierr )
460       ENDIF
461       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
462#endif     
463       IF ( restart_time /= remote )  THEN
464          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
465                 '": restart_time = ', restart_time, '& is not equal to ',     &
466                 'restart_time_remote = ', remote
467          CALL message( 'check_parameters', 'PA0006', 1, 2, 0, 6, 0 )
468       ENDIF
469#if ! defined( __check )
470       IF ( myid == 0 ) THEN
471          CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter,   &
472                         ierr )
473          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter,       &
474                         status, ierr )
475       ENDIF   
476       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
477#endif     
478       IF ( dt_restart /= remote )  THEN
479          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
480                 '": dt_restart = ', dt_restart, '& is not equal to ',         &
481                 'dt_restart_remote = ', remote
482          CALL message( 'check_parameters', 'PA0007', 1, 2, 0, 6, 0 )
483       ENDIF
484
485       simulation_time_since_reference = end_time - coupling_start_time
486#if ! defined( __check )
487       IF  ( myid == 0 ) THEN
488          CALL MPI_SEND( simulation_time_since_reference, 1, MPI_REAL, target_id, &
489                         14, comm_inter, ierr )
490          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter,       &
491                         status, ierr )   
492       ENDIF
493       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
494#endif     
495       IF ( simulation_time_since_reference /= remote )  THEN
496          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
497                 '": simulation_time_since_reference = ',                      &
498                 simulation_time_since_reference, '& is not equal to ',        &
499                 'simulation_time_since_reference_remote = ', remote
500          CALL message( 'check_parameters', 'PA0008', 1, 2, 0, 6, 0 )
501       ENDIF
502
503#if ! defined( __check )
504       IF ( myid == 0 ) THEN
505          CALL MPI_SEND( dx, 1, MPI_REAL, target_id, 15, comm_inter, ierr )
506          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter,       &
507                                                             status, ierr )
508       ENDIF
509       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
510
511#endif
512       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
513
514          IF ( dx < remote ) THEN
515             WRITE( message_string, * ) 'coupling mode "',                     &
516                   TRIM( coupling_mode ),                                      &
517           '": dx in Atmosphere is not equal to or not larger then dx in ocean'
518             CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 )
519          ENDIF
520
521          IF ( (nx_a+1)*dx /= (nx_o+1)*remote )  THEN
522             WRITE( message_string, * ) 'coupling mode "',                     &
523                    TRIM( coupling_mode ),                                     &
524             '": Domain size in x-direction is not equal in ocean and atmosphere'
525             CALL message( 'check_parameters', 'PA0010', 1, 2, 0, 6, 0 )
526          ENDIF
527
528       ENDIF
529
530#if ! defined( __check )
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#endif
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 ) && ! defined ( __check )
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          q_vertical_gradient_level_ind(1) = 0
1546          DO  k = 1, nzt+1
1547             IF ( i < 11 )  THEN
1548                IF ( q_vertical_gradient_level(i) < zu(k)  .AND.               &
1549                     q_vertical_gradient_level(i) >= 0.0_wp )  THEN
1550                   gradient = q_vertical_gradient(i) / 100.0_wp
1551                   q_vertical_gradient_level_ind(i) = k - 1
1552                   i = i + 1
1553                ENDIF
1554             ENDIF
1555             IF ( gradient /= 0.0_wp )  THEN
1556                IF ( k /= 1 )  THEN
1557                   q_init(k) = q_init(k-1) + dzu(k) * gradient
1558                ELSE
1559                   q_init(k) = q_init(k-1) + dzu(k) * gradient
1560                ENDIF
1561             ELSE
1562                q_init(k) = q_init(k-1)
1563             ENDIF
1564!
1565!--          Avoid negative humidities
1566             IF ( q_init(k) < 0.0_wp )  THEN
1567                q_init(k) = 0.0_wp
1568             ENDIF
1569          ENDDO
1570
1571!
1572!--       In case of no given humidity gradients, choose zero gradient
1573!--       conditions
1574          IF ( q_vertical_gradient_level(1) == -1.0_wp )  THEN
1575             q_vertical_gradient_level(1) = 0.0_wp
1576          ENDIF
1577!
1578!--       Store humidity, rain water content and rain drop concentration
1579!--       gradient at the top boundary for possile Neumann boundary condition
1580          bc_q_t_val  = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
1581       ENDIF
1582
1583!
1584!--    If required, compute initial salinity profile using the given salinity
1585!--    gradients
1586       IF ( ocean )  THEN
1587
1588          i = 1
1589          gradient = 0.0_wp
1590
1591          sa_vertical_gradient_level_ind(1) = nzt+1
1592          DO  k = nzt, 0, -1
1593             IF ( i < 11 )  THEN
1594                IF ( sa_vertical_gradient_level(i) > zu(k)  .AND.              &
1595                     sa_vertical_gradient_level(i) <= 0.0_wp )  THEN
1596                   gradient = sa_vertical_gradient(i) / 100.0_wp
1597                   sa_vertical_gradient_level_ind(i) = k + 1
1598                   i = i + 1
1599                ENDIF
1600             ENDIF
1601             IF ( gradient /= 0.0_wp )  THEN
1602                IF ( k /= nzt )  THEN
1603                   sa_init(k) = sa_init(k+1) - dzu(k+1) * gradient
1604                ELSE
1605                   sa_init(k)   = sa_surface - 0.5_wp * dzu(k+1) * gradient
1606                   sa_init(k+1) = sa_surface + 0.5_wp * dzu(k+1) * gradient
1607                ENDIF
1608             ELSE
1609                sa_init(k) = sa_init(k+1)
1610             ENDIF
1611          ENDDO
1612
1613       ENDIF
1614
1615         
1616    ENDIF
1617
1618!
1619!-- Check if the control parameter use_subsidence_tendencies is used correctly
1620    IF ( use_subsidence_tendencies  .AND.  .NOT.  large_scale_subsidence )  THEN
1621       message_string = 'The usage of use_subsidence_tendencies ' //           &
1622                            'requires large_scale_subsidence = .T..'
1623       CALL message( 'check_parameters', 'PA0396', 1, 2, 0, 6, 0 )
1624    ELSEIF ( use_subsidence_tendencies  .AND.  .NOT. large_scale_forcing )  THEN
1625       message_string = 'The usage of use_subsidence_tendencies ' //           &
1626                            'requires large_scale_forcing = .T..'
1627       CALL message( 'check_parameters', 'PA0397', 1, 2, 0, 6, 0 )
1628    ENDIF
1629
1630!
1631!-- Initialize large scale subsidence if required
1632    If ( large_scale_subsidence )  THEN
1633       IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp  .AND.            &
1634                                     .NOT.  large_scale_forcing )  THEN
1635          CALL init_w_subsidence
1636       ENDIF
1637!
1638!--    In case large_scale_forcing is used, profiles for subsidence velocity
1639!--    are read in from file LSF_DATA
1640
1641       IF ( subs_vertical_gradient_level(1) == -9999999.9_wp  .AND.            &
1642            .NOT.  large_scale_forcing )  THEN
1643          message_string = 'There is no default large scale vertical ' //      &
1644                           'velocity profile set. Specify the subsidence ' //  &
1645                           'velocity profile via subs_vertical_gradient ' //   &
1646                           'and subs_vertical_gradient_level.'
1647          CALL message( 'check_parameters', 'PA0380', 1, 2, 0, 6, 0 )
1648       ENDIF
1649    ELSE
1650        IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp )  THEN
1651           message_string = 'Enable usage of large scale subsidence by ' //    &
1652                            'setting large_scale_subsidence = .T..'
1653          CALL message( 'check_parameters', 'PA0381', 1, 2, 0, 6, 0 )
1654        ENDIF
1655    ENDIF   
1656
1657!
1658!-- Compute Coriolis parameter
1659    f  = 2.0_wp * omega * SIN( phi / 180.0_wp * pi )
1660    fs = 2.0_wp * omega * COS( phi / 180.0_wp * pi )
1661
1662!
1663!-- Check and set buoyancy related parameters and switches
1664    IF ( reference_state == 'horizontal_average' )  THEN
1665       CONTINUE
1666    ELSEIF ( reference_state == 'initial_profile' )  THEN
1667       use_initial_profile_as_reference = .TRUE.
1668    ELSEIF ( reference_state == 'single_value' )  THEN
1669       use_single_reference_value = .TRUE.
1670       IF ( pt_reference == 9999999.9_wp )  pt_reference = pt_surface
1671       vpt_reference = pt_reference * ( 1.0_wp + 0.61_wp * q_surface )
1672    ELSE
1673       message_string = 'illegal value for reference_state: "' //              &
1674                        TRIM( reference_state ) // '"'
1675       CALL message( 'check_parameters', 'PA0056', 1, 2, 0, 6, 0 )
1676    ENDIF
1677
1678!
1679!-- Ocean runs always use reference values in the buoyancy term
1680    IF ( ocean )  THEN
1681       reference_state = 'single_value'
1682       use_single_reference_value = .TRUE.
1683    ENDIF
1684
1685!
1686!-- Sign of buoyancy/stability terms
1687    IF ( ocean )  atmos_ocean_sign = -1.0_wp
1688
1689!
1690!-- Ocean version must use flux boundary conditions at the top
1691    IF ( ocean .AND. .NOT. use_top_fluxes )  THEN
1692       message_string = 'use_top_fluxes must be .TRUE. in ocean mode'
1693       CALL message( 'check_parameters', 'PA0042', 1, 2, 0, 6, 0 )
1694    ENDIF
1695
1696!
1697!-- In case of a given slope, compute the relevant quantities
1698    IF ( alpha_surface /= 0.0_wp )  THEN
1699       IF ( ABS( alpha_surface ) > 90.0_wp )  THEN
1700          WRITE( message_string, * ) 'ABS( alpha_surface = ', alpha_surface,   &
1701                                     ' ) must be < 90.0'
1702          CALL message( 'check_parameters', 'PA0043', 1, 2, 0, 6, 0 )
1703       ENDIF
1704       sloping_surface = .TRUE.
1705       cos_alpha_surface = COS( alpha_surface / 180.0_wp * pi )
1706       sin_alpha_surface = SIN( alpha_surface / 180.0_wp * pi )
1707    ENDIF
1708
1709!
1710!-- Check time step and cfl_factor
1711    IF ( dt /= -1.0_wp )  THEN
1712       IF ( dt <= 0.0_wp  .AND.  dt /= -1.0_wp )  THEN
1713          WRITE( message_string, * ) 'dt = ', dt , ' <= 0.0'
1714          CALL message( 'check_parameters', 'PA0044', 1, 2, 0, 6, 0 )
1715       ENDIF
1716       dt_3d = dt
1717       dt_fixed = .TRUE.
1718    ENDIF
1719
1720    IF ( cfl_factor <= 0.0_wp  .OR.  cfl_factor > 1.0_wp )  THEN
1721       IF ( cfl_factor == -1.0_wp )  THEN
1722          IF ( timestep_scheme == 'runge-kutta-2' )  THEN
1723             cfl_factor = 0.8_wp
1724          ELSEIF ( timestep_scheme == 'runge-kutta-3' )  THEN
1725             cfl_factor = 0.9_wp
1726          ELSE
1727             cfl_factor = 0.9_wp
1728          ENDIF
1729       ELSE
1730          WRITE( message_string, * ) 'cfl_factor = ', cfl_factor,              &
1731                 ' out of range & 0.0 < cfl_factor <= 1.0 is required'
1732          CALL message( 'check_parameters', 'PA0045', 1, 2, 0, 6, 0 )
1733       ENDIF
1734    ENDIF
1735
1736!
1737!-- Store simulated time at begin
1738    simulated_time_at_begin = simulated_time
1739
1740!
1741!-- Store reference time for coupled runs and change the coupling flag,
1742!-- if ...
1743    IF ( simulated_time == 0.0_wp )  THEN
1744       IF ( coupling_start_time == 0.0_wp )  THEN
1745          time_since_reference_point = 0.0_wp
1746       ELSEIF ( time_since_reference_point < 0.0_wp )  THEN
1747          run_coupled = .FALSE.
1748       ENDIF
1749    ENDIF
1750
1751!
1752!-- Set wind speed in the Galilei-transformed system
1753    IF ( galilei_transformation )  THEN
1754       IF ( use_ug_for_galilei_tr                    .AND.                     &
1755            ug_vertical_gradient_level(1) == 0.0_wp  .AND.                     &
1756            ug_vertical_gradient(1) == 0.0_wp        .AND.                     & 
1757            vg_vertical_gradient_level(1) == 0.0_wp  .AND.                     &
1758            vg_vertical_gradient(1) == 0.0_wp )  THEN
1759          u_gtrans = ug_surface * 0.6_wp
1760          v_gtrans = vg_surface * 0.6_wp
1761       ELSEIF ( use_ug_for_galilei_tr  .AND.                                   &
1762                ( ug_vertical_gradient_level(1) /= 0.0_wp  .OR.                &
1763                ug_vertical_gradient(1) /= 0.0_wp ) )  THEN
1764          message_string = 'baroclinity (ug) not allowed simultaneously' //    &
1765                           ' with galilei transformation'
1766          CALL message( 'check_parameters', 'PA0046', 1, 2, 0, 6, 0 )
1767       ELSEIF ( use_ug_for_galilei_tr  .AND.                                   &
1768                ( vg_vertical_gradient_level(1) /= 0.0_wp  .OR.                &
1769                vg_vertical_gradient(1) /= 0.0_wp ) )  THEN
1770          message_string = 'baroclinity (vg) not allowed simultaneously' //    &
1771                           ' with galilei transformation'
1772          CALL message( 'check_parameters', 'PA0047', 1, 2, 0, 6, 0 )
1773       ELSE
1774          message_string = 'variable translation speed used for galilei-' //   &
1775             'transformation, which may cause & instabilities in stably ' //   &
1776             'stratified regions'
1777          CALL message( 'check_parameters', 'PA0048', 0, 1, 0, 6, 0 )
1778       ENDIF
1779    ENDIF
1780
1781!
1782!-- In case of using a prandtl-layer, calculated (or prescribed) surface
1783!-- fluxes have to be used in the diffusion-terms
1784    IF ( constant_flux_layer )  use_surface_fluxes = .TRUE.
1785
1786!
1787!-- Check boundary conditions and set internal variables:
1788!-- Attention: the lateral boundary conditions have been already checked in
1789!-- parin
1790!
1791!-- Non-cyclic lateral boundaries require the multigrid method and Piascek-
1792!-- Willimas or Wicker - Skamarock advection scheme. Several schemes
1793!-- and tools do not work with non-cyclic boundary conditions.
1794    IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
1795       IF ( psolver(1:9) /= 'multigrid' )  THEN
1796          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1797                           'psolver = "' // TRIM( psolver ) // '"'
1798          CALL message( 'check_parameters', 'PA0051', 1, 2, 0, 6, 0 )
1799       ENDIF
1800       IF ( momentum_advec /= 'pw-scheme'  .AND.                               &
1801          ( momentum_advec /= 'ws-scheme'  .AND.                               &
1802            momentum_advec /= 'ws-scheme-mono' )                               &
1803          )  THEN
1804
1805          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1806                           'momentum_advec = "' // TRIM( momentum_advec ) // '"'
1807          CALL message( 'check_parameters', 'PA0052', 1, 2, 0, 6, 0 )
1808       ENDIF
1809       IF ( scalar_advec /= 'pw-scheme'  .AND.                                 &
1810          ( scalar_advec /= 'ws-scheme'  .AND.                                 &
1811            scalar_advec /= 'ws-scheme-mono' )                                 &
1812          )  THEN
1813          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1814                           'scalar_advec = "' // TRIM( scalar_advec ) // '"'
1815          CALL message( 'check_parameters', 'PA0053', 1, 2, 0, 6, 0 )
1816       ENDIF
1817       IF ( galilei_transformation )  THEN
1818          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1819                           'galilei_transformation = .T.'
1820          CALL message( 'check_parameters', 'PA0054', 1, 2, 0, 6, 0 )
1821       ENDIF
1822    ENDIF
1823
1824!
1825!-- Bottom boundary condition for the turbulent Kinetic energy
1826    IF ( bc_e_b == 'neumann' )  THEN
1827       ibc_e_b = 1
1828    ELSEIF ( bc_e_b == '(u*)**2+neumann' )  THEN
1829       ibc_e_b = 2
1830       IF ( .NOT. constant_flux_layer )  THEN
1831          bc_e_b = 'neumann'
1832          ibc_e_b = 1
1833          message_string = 'boundary condition bc_e_b changed to "' //         &
1834                           TRIM( bc_e_b ) // '"'
1835          CALL message( 'check_parameters', 'PA0057', 0, 1, 0, 6, 0 )
1836       ENDIF
1837    ELSE
1838       message_string = 'unknown boundary condition: bc_e_b = "' //            &
1839                        TRIM( bc_e_b ) // '"'
1840       CALL message( 'check_parameters', 'PA0058', 1, 2, 0, 6, 0 )
1841    ENDIF
1842
1843!
1844!-- Boundary conditions for perturbation pressure
1845    IF ( bc_p_b == 'dirichlet' )  THEN
1846       ibc_p_b = 0
1847    ELSEIF ( bc_p_b == 'neumann' )  THEN
1848       ibc_p_b = 1
1849    ELSE
1850       message_string = 'unknown boundary condition: bc_p_b = "' //            &
1851                        TRIM( bc_p_b ) // '"'
1852       CALL message( 'check_parameters', 'PA0059', 1, 2, 0, 6, 0 )
1853    ENDIF
1854
1855    IF ( bc_p_t == 'dirichlet' )  THEN
1856       ibc_p_t = 0
1857!-- TO_DO: later set bc_p_t to neumann before, in case of nested domain
1858    ELSEIF ( bc_p_t == 'neumann' .OR. bc_p_t == 'nested' )  THEN
1859       ibc_p_t = 1
1860    ELSE
1861       message_string = 'unknown boundary condition: bc_p_t = "' //            &
1862                        TRIM( bc_p_t ) // '"'
1863       CALL message( 'check_parameters', 'PA0061', 1, 2, 0, 6, 0 )
1864    ENDIF
1865
1866!
1867!-- Boundary conditions for potential temperature
1868    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
1869       ibc_pt_b = 2
1870    ELSE
1871       IF ( bc_pt_b == 'dirichlet' )  THEN
1872          ibc_pt_b = 0
1873       ELSEIF ( bc_pt_b == 'neumann' )  THEN
1874          ibc_pt_b = 1
1875       ELSE
1876          message_string = 'unknown boundary condition: bc_pt_b = "' //        &
1877                           TRIM( bc_pt_b ) // '"'
1878          CALL message( 'check_parameters', 'PA0062', 1, 2, 0, 6, 0 )
1879       ENDIF
1880    ENDIF
1881
1882    IF ( bc_pt_t == 'dirichlet' )  THEN
1883       ibc_pt_t = 0
1884    ELSEIF ( bc_pt_t == 'neumann' )  THEN
1885       ibc_pt_t = 1
1886    ELSEIF ( bc_pt_t == 'initial_gradient' )  THEN
1887       ibc_pt_t = 2
1888    ELSEIF ( bc_pt_t == 'nested' )  THEN
1889       ibc_pt_t = 3
1890    ELSE
1891       message_string = 'unknown boundary condition: bc_pt_t = "' //           &
1892                        TRIM( bc_pt_t ) // '"'
1893       CALL message( 'check_parameters', 'PA0063', 1, 2, 0, 6, 0 )
1894    ENDIF
1895
1896    IF ( surface_heatflux == 9999999.9_wp  )  THEN
1897       constant_heatflux = .FALSE.
1898       IF ( large_scale_forcing  .OR.  land_surface )  THEN
1899          IF ( ibc_pt_b == 0 )  THEN
1900             constant_heatflux = .FALSE.
1901          ELSEIF ( ibc_pt_b == 1 )  THEN
1902             constant_heatflux = .TRUE.
1903             IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.   &
1904                  .NOT.  land_surface )  THEN
1905                surface_heatflux = shf_surf(1)
1906             ELSE
1907                surface_heatflux = 0.0_wp
1908             ENDIF
1909          ENDIF
1910       ENDIF
1911    ELSE
1912        constant_heatflux = .TRUE.
1913        IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.        &
1914             large_scale_forcing  .AND.  .NOT.  land_surface )  THEN
1915           surface_heatflux = shf_surf(1)
1916        ENDIF
1917    ENDIF
1918
1919    IF ( top_heatflux     == 9999999.9_wp )  constant_top_heatflux = .FALSE.
1920
1921    IF ( neutral )  THEN
1922
1923       IF ( surface_heatflux /= 0.0_wp  .AND.                                  &
1924            surface_heatflux /= 9999999.9_wp )  THEN
1925          message_string = 'heatflux must not be set for pure neutral flow'
1926          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
1927       ENDIF
1928
1929       IF ( top_heatflux /= 0.0_wp  .AND.  top_heatflux /= 9999999.9_wp )      &
1930       THEN
1931          message_string = 'heatflux must not be set for pure neutral flow'
1932          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
1933       ENDIF
1934
1935    ENDIF
1936
1937    IF ( top_momentumflux_u /= 9999999.9_wp  .AND.                             &
1938         top_momentumflux_v /= 9999999.9_wp )  THEN
1939       constant_top_momentumflux = .TRUE.
1940    ELSEIF (  .NOT. ( top_momentumflux_u == 9999999.9_wp  .AND.                &
1941           top_momentumflux_v == 9999999.9_wp ) )  THEN
1942       message_string = 'both, top_momentumflux_u AND top_momentumflux_v ' //  &
1943                        'must be set'
1944       CALL message( 'check_parameters', 'PA0064', 1, 2, 0, 6, 0 )
1945    ENDIF
1946
1947!
1948!-- A given surface temperature implies Dirichlet boundary condition for
1949!-- temperature. In this case specification of a constant heat flux is
1950!-- forbidden.
1951    IF ( ibc_pt_b == 0  .AND.  constant_heatflux  .AND.                        &
1952         surface_heatflux /= 0.0_wp )  THEN
1953       message_string = 'boundary_condition: bc_pt_b = "' // TRIM( bc_pt_b ) //&
1954                        '& is not allowed with constant_heatflux = .TRUE.'
1955       CALL message( 'check_parameters', 'PA0065', 1, 2, 0, 6, 0 )
1956    ENDIF
1957    IF ( constant_heatflux  .AND.  pt_surface_initial_change /= 0.0_wp )  THEN
1958       WRITE ( message_string, * )  'constant_heatflux = .TRUE. is not allo',  &
1959               'wed with pt_surface_initial_change (/=0) = ',                  &
1960               pt_surface_initial_change
1961       CALL message( 'check_parameters', 'PA0066', 1, 2, 0, 6, 0 )
1962    ENDIF
1963
1964!
1965!-- A given temperature at the top implies Dirichlet boundary condition for
1966!-- temperature. In this case specification of a constant heat flux is
1967!-- forbidden.
1968    IF ( ibc_pt_t == 0  .AND.  constant_top_heatflux  .AND.                    &
1969         top_heatflux /= 0.0_wp )  THEN
1970       message_string = 'boundary_condition: bc_pt_t = "' // TRIM( bc_pt_t ) //&
1971                        '" is not allowed with constant_top_heatflux = .TRUE.'
1972       CALL message( 'check_parameters', 'PA0067', 1, 2, 0, 6, 0 )
1973    ENDIF
1974
1975!
1976!-- Boundary conditions for salinity
1977    IF ( ocean )  THEN
1978       IF ( bc_sa_t == 'dirichlet' )  THEN
1979          ibc_sa_t = 0
1980       ELSEIF ( bc_sa_t == 'neumann' )  THEN
1981          ibc_sa_t = 1
1982       ELSE
1983          message_string = 'unknown boundary condition: bc_sa_t = "' //        &
1984                           TRIM( bc_sa_t ) // '"'
1985          CALL message( 'check_parameters', 'PA0068', 1, 2, 0, 6, 0 )
1986       ENDIF
1987
1988       IF ( top_salinityflux == 9999999.9_wp )  constant_top_salinityflux = .FALSE.
1989       IF ( ibc_sa_t == 1  .AND.  top_salinityflux == 9999999.9_wp )  THEN
1990          message_string = 'boundary condition: bc_sa_t = "' //                &
1991                           TRIM( bc_sa_t ) // '" requires to set ' //          &
1992                           'top_salinityflux'
1993          CALL message( 'check_parameters', 'PA0069', 1, 2, 0, 6, 0 )
1994       ENDIF
1995
1996!
1997!--    A fixed salinity at the top implies Dirichlet boundary condition for
1998!--    salinity. In this case specification of a constant salinity flux is
1999!--    forbidden.
2000       IF ( ibc_sa_t == 0  .AND.  constant_top_salinityflux  .AND.             &
2001            top_salinityflux /= 0.0_wp )  THEN
2002          message_string = 'boundary condition: bc_sa_t = "' //                &
2003                           TRIM( bc_sa_t ) // '" is not allowed with ' //      &
2004                           'constant_top_salinityflux = .TRUE.'
2005          CALL message( 'check_parameters', 'PA0070', 1, 2, 0, 6, 0 )
2006       ENDIF
2007
2008    ENDIF
2009
2010!
2011!-- In case of humidity or passive scalar, set boundary conditions for total
2012!-- water content / scalar
2013    IF ( humidity  .OR.  passive_scalar ) THEN
2014       IF ( humidity )  THEN
2015          sq = 'q'
2016       ELSE
2017          sq = 's'
2018       ENDIF
2019       IF ( bc_q_b == 'dirichlet' )  THEN
2020          ibc_q_b = 0
2021       ELSEIF ( bc_q_b == 'neumann' )  THEN
2022          ibc_q_b = 1
2023       ELSE
2024          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
2025                           '_b ="' // TRIM( bc_q_b ) // '"'
2026          CALL message( 'check_parameters', 'PA0071', 1, 2, 0, 6, 0 )
2027       ENDIF
2028       IF ( bc_q_t == 'dirichlet' )  THEN
2029          ibc_q_t = 0
2030       ELSEIF ( bc_q_t == 'neumann' )  THEN
2031          ibc_q_t = 1
2032       ELSEIF ( bc_q_t == 'nested' )  THEN
2033          ibc_q_t = 3
2034       ELSE
2035          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
2036                           '_t ="' // TRIM( bc_q_t ) // '"'
2037          CALL message( 'check_parameters', 'PA0072', 1, 2, 0, 6, 0 )
2038       ENDIF
2039
2040       IF ( surface_waterflux == 9999999.9_wp  )  THEN
2041          constant_waterflux = .FALSE.
2042          IF ( large_scale_forcing .OR. land_surface )  THEN
2043             IF ( ibc_q_b == 0 )  THEN
2044                constant_waterflux = .FALSE.
2045             ELSEIF ( ibc_q_b == 1 )  THEN
2046                constant_waterflux = .TRUE.
2047                IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
2048                   surface_waterflux = qsws_surf(1)
2049                ENDIF
2050             ENDIF
2051          ENDIF
2052       ELSE
2053          constant_waterflux = .TRUE.
2054          IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.      &
2055               large_scale_forcing )  THEN
2056             surface_waterflux = qsws_surf(1)
2057          ENDIF
2058       ENDIF
2059
2060!
2061!--    A given surface humidity implies Dirichlet boundary condition for
2062!--    humidity. In this case specification of a constant water flux is
2063!--    forbidden.
2064       IF ( ibc_q_b == 0  .AND.  constant_waterflux )  THEN
2065          message_string = 'boundary condition: bc_' // TRIM( sq ) // '_b ' // &
2066                           '= "' // TRIM( bc_q_b ) // '" is not allowed wi' // &
2067                           'th prescribed surface flux'
2068          CALL message( 'check_parameters', 'PA0073', 1, 2, 0, 6, 0 )
2069       ENDIF
2070       IF ( constant_waterflux  .AND.  q_surface_initial_change /= 0.0_wp )  THEN
2071          WRITE( message_string, * )  'a prescribed surface flux is not allo', &
2072                 'wed with ', sq, '_surface_initial_change (/=0) = ',          &
2073                 q_surface_initial_change
2074          CALL message( 'check_parameters', 'PA0074', 1, 2, 0, 6, 0 )
2075       ENDIF
2076
2077    ENDIF
2078!
2079!-- Boundary conditions for horizontal components of wind speed
2080    IF ( bc_uv_b == 'dirichlet' )  THEN
2081       ibc_uv_b = 0
2082    ELSEIF ( bc_uv_b == 'neumann' )  THEN
2083       ibc_uv_b = 1
2084       IF ( constant_flux_layer )  THEN
2085          message_string = 'boundary condition: bc_uv_b = "' //                &
2086               TRIM( bc_uv_b ) // '" is not allowed with constant_flux_layer'  &
2087               // ' = .TRUE.'
2088          CALL message( 'check_parameters', 'PA0075', 1, 2, 0, 6, 0 )
2089       ENDIF
2090    ELSE
2091       message_string = 'unknown boundary condition: bc_uv_b = "' //           &
2092                        TRIM( bc_uv_b ) // '"'
2093       CALL message( 'check_parameters', 'PA0076', 1, 2, 0, 6, 0 )
2094    ENDIF
2095!
2096!-- In case of coupled simulations u and v at the ground in atmosphere will be
2097!-- assigned with the u and v values of the ocean surface
2098    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
2099       ibc_uv_b = 2
2100    ENDIF
2101
2102    IF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
2103       bc_uv_t = 'neumann'
2104       ibc_uv_t = 1
2105    ELSE
2106       IF ( bc_uv_t == 'dirichlet' .OR. bc_uv_t == 'dirichlet_0' )  THEN
2107          ibc_uv_t = 0
2108          IF ( bc_uv_t == 'dirichlet_0' )  THEN
2109!
2110!--          Velocities for the initial u,v-profiles are set zero at the top
2111!--          in case of dirichlet_0 conditions
2112             u_init(nzt+1)    = 0.0_wp
2113             v_init(nzt+1)    = 0.0_wp
2114          ENDIF
2115       ELSEIF ( bc_uv_t == 'neumann' )  THEN
2116          ibc_uv_t = 1
2117       ELSEIF ( bc_uv_t == 'nested' )  THEN
2118          ibc_uv_t = 3
2119       ELSE
2120          message_string = 'unknown boundary condition: bc_uv_t = "' //        &
2121                           TRIM( bc_uv_t ) // '"'
2122          CALL message( 'check_parameters', 'PA0077', 1, 2, 0, 6, 0 )
2123       ENDIF
2124    ENDIF
2125
2126!
2127!-- Compute and check, respectively, the Rayleigh Damping parameter
2128    IF ( rayleigh_damping_factor == -1.0_wp )  THEN
2129       rayleigh_damping_factor = 0.0_wp
2130    ELSE
2131       IF ( rayleigh_damping_factor < 0.0_wp  .OR.                             &
2132            rayleigh_damping_factor > 1.0_wp )  THEN
2133          WRITE( message_string, * )  'rayleigh_damping_factor = ',            &
2134                              rayleigh_damping_factor, ' out of range [0.0,1.0]'
2135          CALL message( 'check_parameters', 'PA0078', 1, 2, 0, 6, 0 )
2136       ENDIF
2137    ENDIF
2138
2139    IF ( rayleigh_damping_height == -1.0_wp )  THEN
2140       IF (  .NOT.  ocean )  THEN
2141          rayleigh_damping_height = 0.66666666666_wp * zu(nzt)
2142       ELSE
2143          rayleigh_damping_height = 0.66666666666_wp * zu(nzb)
2144       ENDIF
2145    ELSE
2146       IF (  .NOT.  ocean )  THEN
2147          IF ( rayleigh_damping_height < 0.0_wp  .OR.                          &
2148               rayleigh_damping_height > zu(nzt) )  THEN
2149             WRITE( message_string, * )  'rayleigh_damping_height = ',         &
2150                   rayleigh_damping_height, ' out of range [0.0,', zu(nzt), ']'
2151             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
2152          ENDIF
2153       ELSE
2154          IF ( rayleigh_damping_height > 0.0_wp  .OR.                          &
2155               rayleigh_damping_height < zu(nzb) )  THEN
2156             WRITE( message_string, * )  'rayleigh_damping_height = ',         &
2157                   rayleigh_damping_height, ' out of range [0.0,', zu(nzb), ']'
2158             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
2159          ENDIF
2160       ENDIF
2161    ENDIF
2162
2163!
2164!-- Check number of chosen statistic regions. More than 10 regions are not
2165!-- allowed, because so far no more than 10 corresponding output files can
2166!-- be opened (cf. check_open)
2167    IF ( statistic_regions > 9  .OR.  statistic_regions < 0 )  THEN
2168       WRITE ( message_string, * ) 'number of statistic_regions = ',           &
2169                   statistic_regions+1, ' but only 10 regions are allowed'
2170       CALL message( 'check_parameters', 'PA0082', 1, 2, 0, 6, 0 )
2171    ENDIF
2172    IF ( normalizing_region > statistic_regions  .OR.                          &
2173         normalizing_region < 0)  THEN
2174       WRITE ( message_string, * ) 'normalizing_region = ',                    &
2175                normalizing_region, ' must be >= 0 and <= ',statistic_regions, &
2176                ' (value of statistic_regions)'
2177       CALL message( 'check_parameters', 'PA0083', 1, 2, 0, 6, 0 )
2178    ENDIF
2179
2180!
2181!-- Set the default intervals for data output, if necessary
2182!-- NOTE: dt_dosp has already been set in package_parin
2183    IF ( dt_data_output /= 9999999.9_wp )  THEN
2184       IF ( dt_dopr           == 9999999.9_wp )  dt_dopr           = dt_data_output
2185       IF ( dt_dopts          == 9999999.9_wp )  dt_dopts          = dt_data_output
2186       IF ( dt_do2d_xy        == 9999999.9_wp )  dt_do2d_xy        = dt_data_output
2187       IF ( dt_do2d_xz        == 9999999.9_wp )  dt_do2d_xz        = dt_data_output
2188       IF ( dt_do2d_yz        == 9999999.9_wp )  dt_do2d_yz        = dt_data_output
2189       IF ( dt_do3d           == 9999999.9_wp )  dt_do3d           = dt_data_output
2190       IF ( dt_data_output_av == 9999999.9_wp )  dt_data_output_av = dt_data_output
2191       DO  mid = 1, max_masks
2192          IF ( dt_domask(mid) == 9999999.9_wp )  dt_domask(mid)    = dt_data_output
2193       ENDDO
2194    ENDIF
2195
2196!
2197!-- Set the default skip time intervals for data output, if necessary
2198    IF ( skip_time_dopr    == 9999999.9_wp )                                   &
2199                                       skip_time_dopr    = skip_time_data_output
2200    IF ( skip_time_dosp    == 9999999.9_wp )                                   &
2201                                       skip_time_dosp    = skip_time_data_output
2202    IF ( skip_time_do2d_xy == 9999999.9_wp )                                   &
2203                                       skip_time_do2d_xy = skip_time_data_output
2204    IF ( skip_time_do2d_xz == 9999999.9_wp )                                   &
2205                                       skip_time_do2d_xz = skip_time_data_output
2206    IF ( skip_time_do2d_yz == 9999999.9_wp )                                   &
2207                                       skip_time_do2d_yz = skip_time_data_output
2208    IF ( skip_time_do3d    == 9999999.9_wp )                                   &
2209                                       skip_time_do3d    = skip_time_data_output
2210    IF ( skip_time_data_output_av == 9999999.9_wp )                            &
2211                                skip_time_data_output_av = skip_time_data_output
2212    DO  mid = 1, max_masks
2213       IF ( skip_time_domask(mid) == 9999999.9_wp )                            &
2214                                skip_time_domask(mid)    = skip_time_data_output
2215    ENDDO
2216
2217!
2218!-- Check the average intervals (first for 3d-data, then for profiles and
2219!-- spectra)
2220    IF ( averaging_interval > dt_data_output_av )  THEN
2221       WRITE( message_string, * )  'averaging_interval = ',                    &
2222             averaging_interval, ' must be <= dt_data_output = ', dt_data_output
2223       CALL message( 'check_parameters', 'PA0085', 1, 2, 0, 6, 0 )
2224    ENDIF
2225
2226    IF ( averaging_interval_pr == 9999999.9_wp )  THEN
2227       averaging_interval_pr = averaging_interval
2228    ENDIF
2229
2230    IF ( averaging_interval_pr > dt_dopr )  THEN
2231       WRITE( message_string, * )  'averaging_interval_pr = ',                 &
2232             averaging_interval_pr, ' must be <= dt_dopr = ', dt_dopr
2233       CALL message( 'check_parameters', 'PA0086', 1, 2, 0, 6, 0 )
2234    ENDIF
2235
2236    IF ( averaging_interval_sp == 9999999.9_wp )  THEN
2237       averaging_interval_sp = averaging_interval
2238    ENDIF
2239
2240    IF ( averaging_interval_sp > dt_dosp )  THEN
2241       WRITE( message_string, * )  'averaging_interval_sp = ',                 &
2242             averaging_interval_sp, ' must be <= dt_dosp = ', dt_dosp
2243       CALL message( 'check_parameters', 'PA0087', 1, 2, 0, 6, 0 )
2244    ENDIF
2245
2246!
2247!-- Set the default interval for profiles entering the temporal average
2248    IF ( dt_averaging_input_pr == 9999999.9_wp )  THEN
2249       dt_averaging_input_pr = dt_averaging_input
2250    ENDIF
2251
2252!
2253!-- Set the default interval for the output of timeseries to a reasonable
2254!-- value (tries to minimize the number of calls of flow_statistics)
2255    IF ( dt_dots == 9999999.9_wp )  THEN
2256       IF ( averaging_interval_pr == 0.0_wp )  THEN
2257          dt_dots = MIN( dt_run_control, dt_dopr )
2258       ELSE
2259          dt_dots = MIN( dt_run_control, dt_averaging_input_pr )
2260       ENDIF
2261    ENDIF
2262
2263!
2264!-- Check the sample rate for averaging (first for 3d-data, then for profiles)
2265    IF ( dt_averaging_input > averaging_interval )  THEN
2266       WRITE( message_string, * )  'dt_averaging_input = ',                    &
2267                dt_averaging_input, ' must be <= averaging_interval = ',       &
2268                averaging_interval
2269       CALL message( 'check_parameters', 'PA0088', 1, 2, 0, 6, 0 )
2270    ENDIF
2271
2272    IF ( dt_averaging_input_pr > averaging_interval_pr )  THEN
2273       WRITE( message_string, * )  'dt_averaging_input_pr = ',                 &
2274                dt_averaging_input_pr, ' must be <= averaging_interval_pr = ', &
2275                averaging_interval_pr
2276       CALL message( 'check_parameters', 'PA0089', 1, 2, 0, 6, 0 )
2277    ENDIF
2278
2279!
2280!-- Set the default value for the integration interval of precipitation amount
2281    IF ( precipitation )  THEN
2282       IF ( precipitation_amount_interval == 9999999.9_wp )  THEN
2283          precipitation_amount_interval = dt_do2d_xy
2284       ELSE
2285          IF ( precipitation_amount_interval > dt_do2d_xy )  THEN
2286             WRITE( message_string, * )  'precipitation_amount_interval = ',   &
2287                 precipitation_amount_interval, ' must not be larger than ',   &
2288                 'dt_do2d_xy = ', dt_do2d_xy
2289             CALL message( 'check_parameters', 'PA0090', 1, 2, 0, 6, 0 )
2290          ENDIF
2291       ENDIF
2292    ENDIF
2293
2294!
2295!-- Determine the number of output profiles and check whether they are
2296!-- permissible
2297    DO  WHILE ( data_output_pr(dopr_n+1) /= '          ' )
2298
2299       dopr_n = dopr_n + 1
2300       i = dopr_n
2301
2302!
2303!--    Determine internal profile number (for hom, homs)
2304!--    and store height levels
2305       SELECT CASE ( TRIM( data_output_pr(i) ) )
2306
2307          CASE ( 'u', '#u' )
2308             dopr_index(i) = 1
2309             dopr_unit(i)  = 'm/s'
2310             hom(:,2,1,:)  = SPREAD( zu, 2, statistic_regions+1 )
2311             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2312                dopr_initial_index(i) = 5
2313                hom(:,2,5,:)          = SPREAD( zu, 2, statistic_regions+1 )
2314                data_output_pr(i)     = data_output_pr(i)(2:)
2315             ENDIF
2316
2317          CASE ( 'v', '#v' )
2318             dopr_index(i) = 2
2319             dopr_unit(i)  = 'm/s'
2320             hom(:,2,2,:)  = SPREAD( zu, 2, statistic_regions+1 )
2321             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2322                dopr_initial_index(i) = 6
2323                hom(:,2,6,:)          = SPREAD( zu, 2, statistic_regions+1 )
2324                data_output_pr(i)     = data_output_pr(i)(2:)
2325             ENDIF
2326
2327          CASE ( 'w' )
2328             dopr_index(i) = 3
2329             dopr_unit(i)  = 'm/s'
2330             hom(:,2,3,:)  = SPREAD( zw, 2, statistic_regions+1 )
2331
2332          CASE ( 'pt', '#pt' )
2333             IF ( .NOT. cloud_physics ) THEN
2334                dopr_index(i) = 4
2335                dopr_unit(i)  = 'K'
2336                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
2337                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2338                   dopr_initial_index(i) = 7
2339                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
2340                   hom(nzb,2,7,:)        = 0.0_wp    ! because zu(nzb) is negative
2341                   data_output_pr(i)     = data_output_pr(i)(2:)
2342                ENDIF
2343             ELSE
2344                dopr_index(i) = 43
2345                dopr_unit(i)  = 'K'
2346                hom(:,2,43,:)  = SPREAD( zu, 2, statistic_regions+1 )
2347                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2348                   dopr_initial_index(i) = 28
2349                   hom(:,2,28,:)         = SPREAD( zu, 2, statistic_regions+1 )
2350                   hom(nzb,2,28,:)       = 0.0_wp    ! because zu(nzb) is negative
2351                   data_output_pr(i)     = data_output_pr(i)(2:)
2352                ENDIF
2353             ENDIF
2354
2355          CASE ( 'e' )
2356             dopr_index(i)  = 8
2357             dopr_unit(i)   = 'm2/s2'
2358             hom(:,2,8,:)   = SPREAD( zu, 2, statistic_regions+1 )
2359             hom(nzb,2,8,:) = 0.0_wp
2360
2361          CASE ( 'km', '#km' )
2362             dopr_index(i)  = 9
2363             dopr_unit(i)   = 'm2/s'
2364             hom(:,2,9,:)   = SPREAD( zu, 2, statistic_regions+1 )
2365             hom(nzb,2,9,:) = 0.0_wp
2366             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2367                dopr_initial_index(i) = 23
2368                hom(:,2,23,:)         = hom(:,2,9,:)
2369                data_output_pr(i)     = data_output_pr(i)(2:)
2370             ENDIF
2371
2372          CASE ( 'kh', '#kh' )
2373             dopr_index(i)   = 10
2374             dopr_unit(i)    = 'm2/s'
2375             hom(:,2,10,:)   = SPREAD( zu, 2, statistic_regions+1 )
2376             hom(nzb,2,10,:) = 0.0_wp
2377             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2378                dopr_initial_index(i) = 24
2379                hom(:,2,24,:)         = hom(:,2,10,:)
2380                data_output_pr(i)     = data_output_pr(i)(2:)
2381             ENDIF
2382
2383          CASE ( 'l', '#l' )
2384             dopr_index(i)   = 11
2385             dopr_unit(i)    = 'm'
2386             hom(:,2,11,:)   = SPREAD( zu, 2, statistic_regions+1 )
2387             hom(nzb,2,11,:) = 0.0_wp
2388             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2389                dopr_initial_index(i) = 25
2390                hom(:,2,25,:)         = hom(:,2,11,:)
2391                data_output_pr(i)     = data_output_pr(i)(2:)
2392             ENDIF
2393
2394          CASE ( 'w"u"' )
2395             dopr_index(i) = 12
2396             dopr_unit(i)  = 'm2/s2'
2397             hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 )
2398             IF ( constant_flux_layer )  hom(nzb,2,12,:) = zu(1)
2399
2400          CASE ( 'w*u*' )
2401             dopr_index(i) = 13
2402             dopr_unit(i)  = 'm2/s2'
2403             hom(:,2,13,:) = SPREAD( zw, 2, statistic_regions+1 )
2404
2405          CASE ( 'w"v"' )
2406             dopr_index(i) = 14
2407             dopr_unit(i)  = 'm2/s2'
2408             hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 )
2409             IF ( constant_flux_layer )  hom(nzb,2,14,:) = zu(1)
2410
2411          CASE ( 'w*v*' )
2412             dopr_index(i) = 15
2413             dopr_unit(i)  = 'm2/s2'
2414             hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 )
2415
2416          CASE ( 'w"pt"' )
2417             dopr_index(i) = 16
2418             dopr_unit(i)  = 'K m/s'
2419             hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 )
2420
2421          CASE ( 'w*pt*' )
2422             dopr_index(i) = 17
2423             dopr_unit(i)  = 'K m/s'
2424             hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 )
2425
2426          CASE ( 'wpt' )
2427             dopr_index(i) = 18
2428             dopr_unit(i)  = 'K m/s'
2429             hom(:,2,18,:) = SPREAD( zw, 2, statistic_regions+1 )
2430
2431          CASE ( 'wu' )
2432             dopr_index(i) = 19
2433             dopr_unit(i)  = 'm2/s2'
2434             hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 )
2435             IF ( constant_flux_layer )  hom(nzb,2,19,:) = zu(1)
2436
2437          CASE ( 'wv' )
2438             dopr_index(i) = 20
2439             dopr_unit(i)  = 'm2/s2'
2440             hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 )
2441             IF ( constant_flux_layer )  hom(nzb,2,20,:) = zu(1)
2442
2443          CASE ( 'w*pt*BC' )
2444             dopr_index(i) = 21
2445             dopr_unit(i)  = 'K m/s'
2446             hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 )
2447
2448          CASE ( 'wptBC' )
2449             dopr_index(i) = 22
2450             dopr_unit(i)  = 'K m/s'
2451             hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
2452
2453          CASE ( 'sa', '#sa' )
2454             IF ( .NOT. ocean )  THEN
2455                message_string = 'data_output_pr = ' // &
2456                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2457                                 'lemented for ocean = .FALSE.'
2458                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2459             ELSE
2460                dopr_index(i) = 23
2461                dopr_unit(i)  = 'psu'
2462                hom(:,2,23,:) = SPREAD( zu, 2, statistic_regions+1 )
2463                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2464                   dopr_initial_index(i) = 26
2465                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2466                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2467                   data_output_pr(i)     = data_output_pr(i)(2:)
2468                ENDIF
2469             ENDIF
2470
2471          CASE ( 'u*2' )
2472             dopr_index(i) = 30
2473             dopr_unit(i)  = 'm2/s2'
2474             hom(:,2,30,:) = SPREAD( zu, 2, statistic_regions+1 )
2475
2476          CASE ( 'v*2' )
2477             dopr_index(i) = 31
2478             dopr_unit(i)  = 'm2/s2'
2479             hom(:,2,31,:) = SPREAD( zu, 2, statistic_regions+1 )
2480
2481          CASE ( 'w*2' )
2482             dopr_index(i) = 32
2483             dopr_unit(i)  = 'm2/s2'
2484             hom(:,2,32,:) = SPREAD( zw, 2, statistic_regions+1 )
2485
2486          CASE ( 'pt*2' )
2487             dopr_index(i) = 33
2488             dopr_unit(i)  = 'K2'
2489             hom(:,2,33,:) = SPREAD( zu, 2, statistic_regions+1 )
2490
2491          CASE ( 'e*' )
2492             dopr_index(i) = 34
2493             dopr_unit(i)  = 'm2/s2'
2494             hom(:,2,34,:) = SPREAD( zu, 2, statistic_regions+1 )
2495
2496          CASE ( 'w*2pt*' )
2497             dopr_index(i) = 35
2498             dopr_unit(i)  = 'K m2/s2'
2499             hom(:,2,35,:) = SPREAD( zw, 2, statistic_regions+1 )
2500
2501          CASE ( 'w*pt*2' )
2502             dopr_index(i) = 36
2503             dopr_unit(i)  = 'K2 m/s'
2504             hom(:,2,36,:) = SPREAD( zw, 2, statistic_regions+1 )
2505
2506          CASE ( 'w*e*' )
2507             dopr_index(i) = 37
2508             dopr_unit(i)  = 'm3/s3'
2509             hom(:,2,37,:) = SPREAD( zw, 2, statistic_regions+1 )
2510
2511          CASE ( 'w*3' )
2512             dopr_index(i) = 38
2513             dopr_unit(i)  = 'm3/s3'
2514             hom(:,2,38,:) = SPREAD( zw, 2, statistic_regions+1 )
2515
2516          CASE ( 'Sw' )
2517             dopr_index(i) = 39
2518             dopr_unit(i)  = 'none'
2519             hom(:,2,39,:) = SPREAD( zw, 2, statistic_regions+1 )
2520
2521          CASE ( 'p' )
2522             dopr_index(i) = 40
2523             dopr_unit(i)  = 'Pa'
2524             hom(:,2,40,:) = SPREAD( zu, 2, statistic_regions+1 )
2525
2526          CASE ( 'q', '#q' )
2527             IF ( .NOT. humidity )  THEN
2528                message_string = 'data_output_pr = ' // &
2529                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2530                                 'lemented for humidity = .FALSE.'
2531                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2532             ELSE
2533                dopr_index(i) = 41
2534                dopr_unit(i)  = 'kg/kg'
2535                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2536                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2537                   dopr_initial_index(i) = 26
2538                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2539                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2540                   data_output_pr(i)     = data_output_pr(i)(2:)
2541                ENDIF
2542             ENDIF
2543
2544          CASE ( 's', '#s' )
2545             IF ( .NOT. passive_scalar )  THEN
2546                message_string = 'data_output_pr = ' //                        &
2547                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2548                                 'lemented for passive_scalar = .FALSE.'
2549                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2550             ELSE
2551                dopr_index(i) = 41
2552                dopr_unit(i)  = 'kg/m3'
2553                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2554                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2555                   dopr_initial_index(i) = 26
2556                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2557                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2558                   data_output_pr(i)     = data_output_pr(i)(2:)
2559                ENDIF
2560             ENDIF
2561
2562          CASE ( 'qv', '#qv' )
2563             IF ( .NOT. cloud_physics ) THEN
2564                dopr_index(i) = 41
2565                dopr_unit(i)  = 'kg/kg'
2566                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2567                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2568                   dopr_initial_index(i) = 26
2569                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2570                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2571                   data_output_pr(i)     = data_output_pr(i)(2:)
2572                ENDIF
2573             ELSE
2574                dopr_index(i) = 42
2575                dopr_unit(i)  = 'kg/kg'
2576                hom(:,2,42,:) = SPREAD( zu, 2, statistic_regions+1 )
2577                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2578                   dopr_initial_index(i) = 27
2579                   hom(:,2,27,:)         = SPREAD( zu, 2, statistic_regions+1 )
2580                   hom(nzb,2,27,:)       = 0.0_wp   ! because zu(nzb) is negative
2581                   data_output_pr(i)     = data_output_pr(i)(2:)
2582                ENDIF
2583             ENDIF
2584
2585          CASE ( 'lpt', '#lpt' )
2586             IF ( .NOT. cloud_physics ) THEN
2587                message_string = 'data_output_pr = ' //                        &
2588                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2589                                 'lemented for cloud_physics = .FALSE.'
2590                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2591             ELSE
2592                dopr_index(i) = 4
2593                dopr_unit(i)  = 'K'
2594                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
2595                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2596                   dopr_initial_index(i) = 7
2597                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
2598                   hom(nzb,2,7,:)        = 0.0_wp    ! because zu(nzb) is negative
2599                   data_output_pr(i)     = data_output_pr(i)(2:)
2600                ENDIF
2601             ENDIF
2602
2603          CASE ( 'vpt', '#vpt' )
2604             dopr_index(i) = 44
2605             dopr_unit(i)  = 'K'
2606             hom(:,2,44,:) = SPREAD( zu, 2, statistic_regions+1 )
2607             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2608                dopr_initial_index(i) = 29
2609                hom(:,2,29,:)         = SPREAD( zu, 2, statistic_regions+1 )
2610                hom(nzb,2,29,:)       = 0.0_wp    ! because zu(nzb) is negative
2611                data_output_pr(i)     = data_output_pr(i)(2:)
2612             ENDIF
2613
2614          CASE ( 'w"vpt"' )
2615             dopr_index(i) = 45
2616             dopr_unit(i)  = 'K m/s'
2617             hom(:,2,45,:) = SPREAD( zw, 2, statistic_regions+1 )
2618
2619          CASE ( 'w*vpt*' )
2620             dopr_index(i) = 46
2621             dopr_unit(i)  = 'K m/s'
2622             hom(:,2,46,:) = SPREAD( zw, 2, statistic_regions+1 )
2623
2624          CASE ( 'wvpt' )
2625             dopr_index(i) = 47
2626             dopr_unit(i)  = 'K m/s'
2627             hom(:,2,47,:) = SPREAD( zw, 2, statistic_regions+1 )
2628
2629          CASE ( 'w"q"' )
2630             IF (  .NOT.  humidity )  THEN
2631                message_string = 'data_output_pr = ' //                        &
2632                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2633                                 'lemented for humidity = .FALSE.'
2634                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2635             ELSE
2636                dopr_index(i) = 48
2637                dopr_unit(i)  = 'kg/kg m/s'
2638                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2639             ENDIF
2640
2641          CASE ( 'w*q*' )
2642             IF (  .NOT.  humidity )  THEN
2643                message_string = 'data_output_pr = ' //                        &
2644                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2645                                 'lemented for humidity = .FALSE.'
2646                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2647             ELSE
2648                dopr_index(i) = 49
2649                dopr_unit(i)  = 'kg/kg m/s'
2650                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2651             ENDIF
2652
2653          CASE ( 'wq' )
2654             IF (  .NOT.  humidity )  THEN
2655                message_string = 'data_output_pr = ' //                        &
2656                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2657                                 'lemented for humidity = .FALSE.'
2658                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2659             ELSE
2660                dopr_index(i) = 50
2661                dopr_unit(i)  = 'kg/kg m/s'
2662                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2663             ENDIF
2664
2665          CASE ( 'w"s"' )
2666             IF (  .NOT.  passive_scalar )  THEN
2667                message_string = 'data_output_pr = ' //                        &
2668                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2669                                 'lemented for passive_scalar = .FALSE.'
2670                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2671             ELSE
2672                dopr_index(i) = 48
2673                dopr_unit(i)  = 'kg/m3 m/s'
2674                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2675             ENDIF
2676
2677          CASE ( 'w*s*' )
2678             IF (  .NOT.  passive_scalar )  THEN
2679                message_string = 'data_output_pr = ' //                        &
2680                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2681                                 'lemented for passive_scalar = .FALSE.'
2682                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2683             ELSE
2684                dopr_index(i) = 49
2685                dopr_unit(i)  = 'kg/m3 m/s'
2686                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2687             ENDIF
2688
2689          CASE ( 'ws' )
2690             IF (  .NOT.  passive_scalar )  THEN
2691                message_string = 'data_output_pr = ' //                        &
2692                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2693                                 'lemented for passive_scalar = .FALSE.'
2694                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2695             ELSE
2696                dopr_index(i) = 50
2697                dopr_unit(i)  = 'kg/m3 m/s'
2698                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2699             ENDIF
2700
2701          CASE ( 'w"qv"' )
2702             IF ( humidity  .AND.  .NOT.  cloud_physics )  THEN
2703                dopr_index(i) = 48
2704                dopr_unit(i)  = 'kg/kg m/s'
2705                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2706             ELSEIF ( humidity  .AND.  cloud_physics )  THEN
2707                dopr_index(i) = 51
2708                dopr_unit(i)  = 'kg/kg m/s'
2709                hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
2710             ELSE
2711                message_string = 'data_output_pr = ' //                        &
2712                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2713                                 'lemented for cloud_physics = .FALSE. an&' // &
2714                                 'd humidity = .FALSE.'
2715                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2716             ENDIF
2717
2718          CASE ( 'w*qv*' )
2719             IF ( humidity  .AND.  .NOT. cloud_physics )                       &
2720             THEN
2721                dopr_index(i) = 49
2722                dopr_unit(i)  = 'kg/kg m/s'
2723                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2724             ELSEIF( humidity .AND. cloud_physics ) THEN
2725                dopr_index(i) = 52
2726                dopr_unit(i)  = 'kg/kg m/s'
2727                hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
2728             ELSE
2729                message_string = 'data_output_pr = ' //                        &
2730                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2731                                 'lemented for cloud_physics = .FALSE. an&' // &
2732                                 'd humidity = .FALSE.'
2733                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2734             ENDIF
2735
2736          CASE ( 'wqv' )
2737             IF ( humidity  .AND.  .NOT.  cloud_physics )  THEN
2738                dopr_index(i) = 50
2739                dopr_unit(i)  = 'kg/kg m/s'
2740                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2741             ELSEIF ( humidity  .AND.  cloud_physics )  THEN
2742                dopr_index(i) = 53
2743                dopr_unit(i)  = 'kg/kg m/s'
2744                hom(:,2,53,:) = SPREAD( zw, 2, statistic_regions+1 )
2745             ELSE
2746                message_string = 'data_output_pr = ' //                        &
2747                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2748                                 'lemented for cloud_physics = .FALSE. an&' // &
2749                                 'd humidity = .FALSE.'
2750                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2751             ENDIF
2752
2753          CASE ( 'ql' )
2754             IF (  .NOT.  cloud_physics  .AND.  .NOT.  cloud_droplets )  THEN
2755                message_string = 'data_output_pr = ' //                        &
2756                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2757                                 'lemented for cloud_physics = .FALSE. or'  // &
2758                                 '&cloud_droplets = .FALSE.'
2759                CALL message( 'check_parameters', 'PA0096', 1, 2, 0, 6, 0 )
2760             ELSE
2761                dopr_index(i) = 54
2762                dopr_unit(i)  = 'kg/kg'
2763                hom(:,2,54,:)  = SPREAD( zu, 2, statistic_regions+1 )
2764             ENDIF
2765
2766          CASE ( 'w*u*u*:dz' )
2767             dopr_index(i) = 55
2768             dopr_unit(i)  = 'm2/s3'
2769             hom(:,2,55,:) = SPREAD( zu, 2, statistic_regions+1 )
2770
2771          CASE ( 'w*p*:dz' )
2772             dopr_index(i) = 56
2773             dopr_unit(i)  = 'm2/s3'
2774             hom(:,2,56,:) = SPREAD( zw, 2, statistic_regions+1 )
2775
2776          CASE ( 'w"e:dz' )
2777             dopr_index(i) = 57
2778             dopr_unit(i)  = 'm2/s3'
2779             hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 )
2780
2781
2782          CASE ( 'u"pt"' )
2783             dopr_index(i) = 58
2784             dopr_unit(i)  = 'K m/s'
2785             hom(:,2,58,:) = SPREAD( zu, 2, statistic_regions+1 )
2786
2787          CASE ( 'u*pt*' )
2788             dopr_index(i) = 59
2789             dopr_unit(i)  = 'K m/s'
2790             hom(:,2,59,:) = SPREAD( zu, 2, statistic_regions+1 )
2791
2792          CASE ( 'upt_t' )
2793             dopr_index(i) = 60
2794             dopr_unit(i)  = 'K m/s'
2795             hom(:,2,60,:) = SPREAD( zu, 2, statistic_regions+1 )
2796
2797          CASE ( 'v"pt"' )
2798             dopr_index(i) = 61
2799             dopr_unit(i)  = 'K m/s'
2800             hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
2801             
2802          CASE ( 'v*pt*' )
2803             dopr_index(i) = 62
2804             dopr_unit(i)  = 'K m/s'
2805             hom(:,2,62,:) = SPREAD( zu, 2, statistic_regions+1 )
2806
2807          CASE ( 'vpt_t' )
2808             dopr_index(i) = 63
2809             dopr_unit(i)  = 'K m/s'
2810             hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
2811
2812          CASE ( 'rho' )
2813             IF (  .NOT.  ocean ) THEN
2814                message_string = 'data_output_pr = ' //                        &
2815                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2816                                 'lemented for ocean = .FALSE.'
2817                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2818             ELSE
2819                dopr_index(i) = 64
2820                dopr_unit(i)  = 'kg/m3'
2821                hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
2822                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2823                   dopr_initial_index(i) = 77
2824                   hom(:,2,77,:)         = SPREAD( zu, 2, statistic_regions+1 )
2825                   hom(nzb,2,77,:)       = 0.0_wp    ! because zu(nzb) is negative
2826                   data_output_pr(i)     = data_output_pr(i)(2:)
2827                ENDIF
2828             ENDIF
2829
2830          CASE ( 'w"sa"' )
2831             IF (  .NOT.  ocean ) THEN
2832                message_string = 'data_output_pr = ' //                        &
2833                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2834                                 'lemented for ocean = .FALSE.'
2835                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2836             ELSE
2837                dopr_index(i) = 65
2838                dopr_unit(i)  = 'psu m/s'
2839                hom(:,2,65,:) = SPREAD( zw, 2, statistic_regions+1 )
2840             ENDIF
2841
2842          CASE ( 'w*sa*' )
2843             IF (  .NOT. ocean  ) THEN
2844                message_string = 'data_output_pr = ' //                        &
2845                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2846                                 'lemented for ocean = .FALSE.'
2847                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2848             ELSE
2849                dopr_index(i) = 66
2850                dopr_unit(i)  = 'psu m/s'
2851                hom(:,2,66,:) = SPREAD( zw, 2, statistic_regions+1 )
2852             ENDIF
2853
2854          CASE ( 'wsa' )
2855             IF (  .NOT.  ocean ) THEN
2856                message_string = 'data_output_pr = ' //                        &
2857                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2858                                 'lemented for ocean = .FALSE.'
2859                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2860             ELSE
2861                dopr_index(i) = 67
2862                dopr_unit(i)  = 'psu m/s'
2863                hom(:,2,67,:) = SPREAD( zw, 2, statistic_regions+1 )
2864             ENDIF
2865
2866          CASE ( 'w*p*' )
2867             dopr_index(i) = 68
2868             dopr_unit(i)  = 'm3/s3'
2869             hom(:,2,68,:) = SPREAD( zu, 2, statistic_regions+1 )
2870
2871          CASE ( 'w"e' )
2872             dopr_index(i) = 69
2873             dopr_unit(i)  = 'm3/s3'
2874             hom(:,2,69,:) = SPREAD( zu, 2, statistic_regions+1 )
2875
2876          CASE ( 'q*2' )
2877             IF (  .NOT.  humidity )  THEN
2878                message_string = 'data_output_pr = ' //                        &
2879                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2880                                 'lemented for humidity = .FALSE.'
2881                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2882             ELSE
2883                dopr_index(i) = 70
2884                dopr_unit(i)  = 'kg2/kg2'
2885                hom(:,2,70,:) = SPREAD( zu, 2, statistic_regions+1 )
2886             ENDIF
2887
2888          CASE ( 'prho' )
2889             IF (  .NOT.  ocean ) THEN
2890                message_string = 'data_output_pr = ' //                        &
2891                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2892                                 'lemented for ocean = .FALSE.'
2893                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2894             ELSE
2895                dopr_index(i) = 71
2896                dopr_unit(i)  = 'kg/m3'
2897                hom(:,2,71,:) = SPREAD( zu, 2, statistic_regions+1 )
2898             ENDIF
2899
2900          CASE ( 'hyp' )
2901             dopr_index(i) = 72
2902             dopr_unit(i)  = 'dbar'
2903             hom(:,2,72,:) = SPREAD( zu, 2, statistic_regions+1 )
2904
2905          CASE ( 'nr' )
2906             IF (  .NOT.  cloud_physics )  THEN
2907                message_string = 'data_output_pr = ' //                        &
2908                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2909                                 'lemented for cloud_physics = .FALSE.'
2910                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2911             ELSEIF ( icloud_scheme /= 0 )  THEN
2912                message_string = 'data_output_pr = ' //                        &
2913                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2914                                 'lemented for cloud_scheme /= seifert_beheng'
2915                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2916             ELSEIF (  .NOT.  precipitation )  THEN
2917                message_string = 'data_output_pr = ' //                        &
2918                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2919                                 'lemented for precipitation = .FALSE.'
2920                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
2921             ELSE
2922                dopr_index(i) = 73
2923                dopr_unit(i)  = '1/m3'
2924                hom(:,2,73,:)  = SPREAD( zu, 2, statistic_regions+1 )
2925             ENDIF
2926
2927          CASE ( 'qr' )
2928             IF (  .NOT.  cloud_physics )  THEN
2929                message_string = 'data_output_pr = ' //                        &
2930                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2931                                 'lemented for cloud_physics = .FALSE.'
2932                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2933             ELSEIF ( icloud_scheme /= 0 )  THEN
2934                message_string = 'data_output_pr = ' //                        &
2935                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2936                                 'lemented for cloud_scheme /= seifert_beheng'
2937                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2938             ELSEIF (  .NOT.  precipitation )  THEN
2939                message_string = 'data_output_pr = ' //                        &
2940                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2941                                 'lemented for precipitation = .FALSE.'
2942                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
2943             ELSE
2944                dopr_index(i) = 74
2945                dopr_unit(i)  = 'kg/kg'
2946                hom(:,2,74,:)  = SPREAD( zu, 2, statistic_regions+1 )
2947             ENDIF
2948
2949          CASE ( 'qc' )
2950             IF (  .NOT.  cloud_physics )  THEN
2951                message_string = 'data_output_pr = ' //                        &
2952                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2953                                 'lemented for cloud_physics = .FALSE.'
2954                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2955             ELSEIF ( icloud_scheme /= 0 )  THEN
2956                message_string = 'data_output_pr = ' //                        &
2957                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2958                                 'lemented for cloud_scheme /= seifert_beheng'
2959                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2960             ELSE
2961                dopr_index(i) = 75
2962                dopr_unit(i)  = 'kg/kg'
2963                hom(:,2,75,:)  = SPREAD( zu, 2, statistic_regions+1 )
2964             ENDIF
2965
2966          CASE ( 'prr' )
2967             IF (  .NOT.  cloud_physics )  THEN
2968                message_string = 'data_output_pr = ' //                        &
2969                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2970                                 'lemented for cloud_physics = .FALSE.'
2971                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2972             ELSEIF ( icloud_scheme /= 0 )  THEN
2973                message_string = 'data_output_pr = ' //                        &
2974                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2975                                 'lemented for cloud_scheme /= seifert_beheng'
2976                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2977             ELSEIF (  .NOT.  precipitation )  THEN
2978                message_string = 'data_output_pr = ' //                        &
2979                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2980                                 'lemented for precipitation = .FALSE.'
2981                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
2982
2983             ELSE
2984                dopr_index(i) = 76
2985                dopr_unit(i)  = 'kg/kg m/s'
2986                hom(:,2,76,:)  = SPREAD( zu, 2, statistic_regions+1 )
2987             ENDIF
2988
2989          CASE ( 'ug' )
2990             dopr_index(i) = 78
2991             dopr_unit(i)  = 'm/s'
2992             hom(:,2,78,:) = SPREAD( zu, 2, statistic_regions+1 )
2993
2994          CASE ( 'vg' )
2995             dopr_index(i) = 79
2996             dopr_unit(i)  = 'm/s'
2997             hom(:,2,79,:) = SPREAD( zu, 2, statistic_regions+1 )
2998
2999          CASE ( 'w_subs' )
3000             IF (  .NOT.  large_scale_subsidence )  THEN
3001                message_string = 'data_output_pr = ' //                        &
3002                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3003                                 'lemented for large_scale_subsidence = .FALSE.'
3004                CALL message( 'check_parameters', 'PA0382', 1, 2, 0, 6, 0 )
3005             ELSE
3006                dopr_index(i) = 80
3007                dopr_unit(i)  = 'm/s'
3008                hom(:,2,80,:) = SPREAD( zu, 2, statistic_regions+1 )
3009             ENDIF
3010
3011          CASE ( 'td_lsa_lpt' )
3012             IF (  .NOT.  large_scale_forcing )  THEN
3013                message_string = 'data_output_pr = ' //                        &
3014                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3015                                 'lemented for large_scale_forcing = .FALSE.'
3016                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
3017             ELSE
3018                dopr_index(i) = 81
3019                dopr_unit(i)  = 'K/s'
3020                hom(:,2,81,:) = SPREAD( zu, 2, statistic_regions+1 )
3021             ENDIF
3022
3023          CASE ( 'td_lsa_q' )
3024             IF (  .NOT.  large_scale_forcing )  THEN
3025                message_string = 'data_output_pr = ' //                        &
3026                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3027                                 'lemented for large_scale_forcing = .FALSE.'
3028                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
3029             ELSE
3030                dopr_index(i) = 82
3031                dopr_unit(i)  = 'kg/kgs'
3032                hom(:,2,82,:) = SPREAD( zu, 2, statistic_regions+1 )
3033             ENDIF
3034
3035          CASE ( 'td_sub_lpt' )
3036             IF (  .NOT.  large_scale_forcing )  THEN
3037                message_string = 'data_output_pr = ' //                        &
3038                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3039                                 'lemented for large_scale_forcing = .FALSE.'
3040                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
3041             ELSE
3042                dopr_index(i) = 83
3043                dopr_unit(i)  = 'K/s'
3044                hom(:,2,83,:) = SPREAD( zu, 2, statistic_regions+1 )
3045             ENDIF
3046
3047          CASE ( 'td_sub_q' )
3048             IF (  .NOT.  large_scale_forcing )  THEN
3049                message_string = 'data_output_pr = ' //                        &
3050                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3051                                 'lemented for large_scale_forcing = .FALSE.'
3052                CALL message( 'check_parameters', 'PA0393', 1, 2, 0, 6, 0 )
3053             ELSE
3054                dopr_index(i) = 84
3055                dopr_unit(i)  = 'kg/kgs'
3056                hom(:,2,84,:) = SPREAD( zu, 2, statistic_regions+1 )
3057             ENDIF
3058
3059          CASE ( 'td_nud_lpt' )
3060             IF (  .NOT.  nudging )  THEN
3061                message_string = 'data_output_pr = ' //                        &
3062                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3063                                 'lemented for nudging = .FALSE.'
3064                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
3065             ELSE
3066                dopr_index(i) = 85
3067                dopr_unit(i)  = 'K/s'
3068                hom(:,2,85,:) = SPREAD( zu, 2, statistic_regions+1 )
3069             ENDIF
3070
3071          CASE ( 'td_nud_q' )
3072             IF (  .NOT.  nudging )  THEN
3073                message_string = 'data_output_pr = ' //                        &
3074                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3075                                 'lemented for nudging = .FALSE.'
3076                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
3077             ELSE
3078                dopr_index(i) = 86
3079                dopr_unit(i)  = 'kg/kgs'
3080                hom(:,2,86,:) = SPREAD( zu, 2, statistic_regions+1 )
3081             ENDIF
3082
3083          CASE ( 'td_nud_u' )
3084             IF (  .NOT.  nudging )  THEN
3085                message_string = 'data_output_pr = ' //                        &
3086                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3087                                 'lemented for nudging = .FALSE.'
3088                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
3089             ELSE
3090                dopr_index(i) = 87
3091                dopr_unit(i)  = 'm/s2'
3092                hom(:,2,87,:) = SPREAD( zu, 2, statistic_regions+1 )
3093             ENDIF
3094
3095          CASE ( 'td_nud_v' )
3096             IF (  .NOT.  nudging )  THEN
3097                message_string = 'data_output_pr = ' //                        &
3098                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3099                                 'lemented for nudging = .FALSE.'
3100                CALL message( 'check_parameters', 'PA0394', 1, 2, 0, 6, 0 )
3101             ELSE
3102                dopr_index(i) = 88
3103                dopr_unit(i)  = 'm/s2'
3104                hom(:,2,88,:) = SPREAD( zu, 2, statistic_regions+1 )
3105             ENDIF
3106
3107          CASE ( 't_soil', '#t_soil' )
3108             IF (  .NOT.  land_surface )  THEN
3109                message_string = 'data_output_pr = ' //                        &
3110                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3111                                 'lemented for land_surface = .FALSE.'
3112                CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
3113             ELSE
3114                dopr_index(i) = 89
3115                dopr_unit(i)  = 'K'
3116                hom(0:nzs-1,2,89,:)  = SPREAD( - zs, 2, statistic_regions+1 )
3117                IF ( data_output_pr(i)(1:1) == '#' )  THEN
3118                   dopr_initial_index(i) = 90
3119                   hom(0:nzs-1,2,90,:)   = SPREAD( - zs, 2, statistic_regions+1 )
3120                   data_output_pr(i)     = data_output_pr(i)(2:)
3121                ENDIF
3122             ENDIF
3123
3124          CASE ( 'm_soil', '#m_soil' )
3125             IF (  .NOT.  land_surface )  THEN
3126                message_string = 'data_output_pr = ' //                        &
3127                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3128                                 'lemented for land_surface = .FALSE.'
3129                CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
3130             ELSE
3131                dopr_index(i) = 91
3132                dopr_unit(i)  = 'm3/m3'
3133                hom(0:nzs-1,2,91,:)  = SPREAD( - zs, 2, statistic_regions+1 )
3134                IF ( data_output_pr(i)(1:1) == '#' )  THEN
3135                   dopr_initial_index(i) = 92
3136                   hom(0:nzs-1,2,92,:)   = SPREAD( - zs, 2, statistic_regions+1 )
3137                   data_output_pr(i)     = data_output_pr(i)(2:)
3138                ENDIF
3139             ENDIF
3140
3141          CASE ( 'rad_net' )
3142             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
3143             THEN
3144                message_string = 'data_output_pr = ' //                        &
3145                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
3146                                 'lable for radiation = .FALSE. or ' //        &
3147                                 'radiation_scheme = "constant"'
3148                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
3149             ELSE
3150                dopr_index(i) = 101
3151                dopr_unit(i)  = 'W/m2'
3152                hom(:,2,101,:)  = SPREAD( zw, 2, statistic_regions+1 )
3153             ENDIF
3154
3155          CASE ( 'rad_lw_in' )
3156             IF ( (  .NOT.  radiation)  .OR.  radiation_scheme == 'constant' ) &
3157             THEN
3158                message_string = 'data_output_pr = ' //                        &
3159                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
3160                                 'lable for radiation = .FALSE. or ' //        &
3161                                 'radiation_scheme = "constant"'
3162                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
3163             ELSE
3164                dopr_index(i) = 102
3165                dopr_unit(i)  = 'W/m2'
3166                hom(:,2,102,:)  = SPREAD( zw, 2, statistic_regions+1 )
3167             ENDIF
3168
3169          CASE ( 'rad_lw_out' )
3170             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
3171             THEN
3172                message_string = 'data_output_pr = ' //                        &
3173                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
3174                                 'lable for radiation = .FALSE. or ' //        &
3175                                 'radiation_scheme = "constant"'
3176                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
3177             ELSE
3178                dopr_index(i) = 103
3179                dopr_unit(i)  = 'W/m2'
3180                hom(:,2,103,:)  = SPREAD( zw, 2, statistic_regions+1 )
3181             ENDIF
3182
3183          CASE ( 'rad_sw_in' )
3184             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
3185             THEN
3186                message_string = 'data_output_pr = ' //                        &
3187                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
3188                                 'lable for radiation = .FALSE. or ' //        &
3189                                 'radiation_scheme = "constant"'
3190                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
3191             ELSE
3192                dopr_index(i) = 104
3193                dopr_unit(i)  = 'W/m2'
3194                hom(:,2,104,:)  = SPREAD( zw, 2, statistic_regions+1 )
3195             ENDIF
3196
3197          CASE ( 'rad_sw_out')
3198             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' ) &
3199             THEN
3200                message_string = 'data_output_pr = ' //                        &
3201                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
3202                                 'lable for radiation = .FALSE. or ' //        &
3203                                 'radiation_scheme = "constant"'
3204                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
3205             ELSE
3206                dopr_index(i) = 105
3207                dopr_unit(i)  = 'W/m2'
3208                hom(:,2,105,:)  = SPREAD( zw, 2, statistic_regions+1 )
3209             ENDIF
3210
3211          CASE ( 'rad_lw_cs_hr' )
3212             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
3213             THEN
3214                message_string = 'data_output_pr = ' //                        &
3215                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
3216                                 'lable for radiation = .FALSE. or ' //        &
3217                                 'radiation_scheme /= "rrtmg"'
3218                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
3219             ELSE
3220                dopr_index(i) = 106
3221                dopr_unit(i)  = 'K/h'
3222                hom(:,2,106,:)  = SPREAD( zu, 2, statistic_regions+1 )
3223             ENDIF
3224
3225          CASE ( 'rad_lw_hr' )
3226             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
3227             THEN
3228                message_string = 'data_output_pr = ' //                        &
3229                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
3230                                 'lable for radiation = .FALSE. or ' //        &
3231                                 'radiation_scheme /= "rrtmg"'
3232                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
3233             ELSE
3234                dopr_index(i) = 107
3235                dopr_unit(i)  = 'K/h'
3236                hom(:,2,107,:)  = SPREAD( zu, 2, statistic_regions+1 )
3237             ENDIF
3238
3239          CASE ( 'rad_sw_cs_hr' )
3240             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
3241             THEN
3242                message_string = 'data_output_pr = ' //                        &
3243                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
3244                                 'lable for radiation = .FALSE. or ' //        &
3245                                 'radiation_scheme /= "rrtmg"'
3246                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
3247             ELSE
3248                dopr_index(i) = 108
3249                dopr_unit(i)  = 'K/h'
3250                hom(:,2,108,:)  = SPREAD( zu, 2, statistic_regions+1 )
3251             ENDIF
3252
3253          CASE ( 'rad_sw_hr' )
3254             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
3255             THEN
3256                message_string = 'data_output_pr = ' //                        &
3257                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
3258                                 'lable for radiation = .FALSE. or ' //        &
3259                                 'radiation_scheme /= "rrtmg"'
3260                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
3261             ELSE
3262                dopr_index(i) = 109
3263                dopr_unit(i)  = 'K/h'
3264                hom(:,2,109,:)  = SPREAD( zu, 2, statistic_regions+1 )
3265             ENDIF
3266
3267          CASE DEFAULT
3268
3269             CALL user_check_data_output_pr( data_output_pr(i), i, unit )
3270
3271             IF ( unit == 'illegal' )  THEN
3272                IF ( data_output_pr_user(1) /= ' ' )  THEN
3273                   message_string = 'illegal value for data_output_pr or ' //  &
3274                                    'data_output_pr_user = "' //               &
3275                                    TRIM( data_output_pr(i) ) // '"'
3276                   CALL message( 'check_parameters', 'PA0097', 1, 2, 0, 6, 0 )
3277                ELSE
3278                   message_string = 'illegal value for data_output_pr = "' //  &
3279                                    TRIM( data_output_pr(i) ) // '"'
3280                   CALL message( 'check_parameters', 'PA0098', 1, 2, 0, 6, 0 )
3281                ENDIF
3282             ENDIF
3283
3284       END SELECT
3285
3286    ENDDO
3287
3288
3289!
3290!-- Append user-defined data output variables to the standard data output
3291    IF ( data_output_user(1) /= ' ' )  THEN
3292       i = 1
3293       DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 100 )
3294          i = i + 1
3295       ENDDO
3296       j = 1
3297       DO  WHILE ( data_output_user(j) /= ' '  .AND.  j <= 100 )
3298          IF ( i > 100 )  THEN
3299             message_string = 'number of output quantitities given by data' // &
3300                '_output and data_output_user exceeds the limit of 100'
3301             CALL message( 'check_parameters', 'PA0102', 1, 2, 0, 6, 0 )
3302          ENDIF
3303          data_output(i) = data_output_user(j)
3304          i = i + 1
3305          j = j + 1
3306       ENDDO
3307    ENDIF
3308
3309!
3310!-- Check and set steering parameters for 2d/3d data output and averaging
3311    i   = 1
3312    DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 100 )
3313!
3314!--    Check for data averaging
3315       ilen = LEN_TRIM( data_output(i) )
3316       j = 0                                                 ! no data averaging
3317       IF ( ilen > 3 )  THEN
3318          IF ( data_output(i)(ilen-2:ilen) == '_av' )  THEN
3319             j = 1                                           ! data averaging
3320             data_output(i) = data_output(i)(1:ilen-3)
3321          ENDIF
3322       ENDIF
3323!
3324!--    Check for cross section or volume data
3325       ilen = LEN_TRIM( data_output(i) )
3326       k = 0                                                   ! 3d data
3327       var = data_output(i)(1:ilen)
3328       IF ( ilen > 3 )  THEN
3329          IF ( data_output(i)(ilen-2:ilen) == '_xy'  .OR.                      &
3330               data_output(i)(ilen-2:ilen) == '_xz'  .OR.                      &
3331               data_output(i)(ilen-2:ilen) == '_yz' )  THEN
3332             k = 1                                             ! 2d data
3333             var = data_output(i)(1:ilen-3)
3334          ENDIF
3335       ENDIF
3336
3337!
3338!--    Check for allowed value and set units
3339       SELECT CASE ( TRIM( var ) )
3340
3341          CASE ( 'e' )
3342             IF ( constant_diffusion )  THEN
3343                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3344                                 'res constant_diffusion = .FALSE.'
3345                CALL message( 'check_parameters', 'PA0103', 1, 2, 0, 6, 0 )
3346             ENDIF
3347             unit = 'm2/s2'
3348
3349          CASE ( 'lpt' )
3350             IF (  .NOT.  cloud_physics )  THEN
3351                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3352                         'res cloud_physics = .TRUE.'
3353                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3354             ENDIF
3355             unit = 'K'
3356
3357          CASE ( 'm_soil' )
3358             IF (  .NOT.  land_surface )  THEN
3359                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3360                         'land_surface = .TRUE.'
3361                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3362             ENDIF
3363             unit = 'm3/m3'
3364
3365          CASE ( 'nr' )
3366             IF (  .NOT.  cloud_physics )  THEN
3367                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3368                         'res cloud_physics = .TRUE.'
3369                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3370             ELSEIF ( icloud_scheme /= 0 )  THEN
3371                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3372                         'res cloud_scheme = seifert_beheng'
3373                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
3374             ENDIF
3375             unit = '1/m3'
3376
3377          CASE ( 'pc', 'pr' )
3378             IF (  .NOT.  particle_advection )  THEN
3379                message_string = 'output of "' // TRIM( var ) // '" requir' // &
3380                   'es a "particles_par"-NAMELIST in the parameter file (PARIN)'
3381                CALL message( 'check_parameters', 'PA0104', 1, 2, 0, 6, 0 )
3382             ENDIF
3383             IF ( TRIM( var ) == 'pc' )  unit = 'number'
3384             IF ( TRIM( var ) == 'pr' )  unit = 'm'
3385
3386          CASE ( 'prr' )
3387             IF (  .NOT.  cloud_physics )  THEN
3388                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3389                         'res cloud_physics = .TRUE.'
3390                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3391             ELSEIF ( icloud_scheme /= 0 )  THEN
3392                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3393                         'res cloud_scheme = seifert_beheng'
3394                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
3395             ELSEIF (  .NOT.  precipitation )  THEN
3396                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3397                                 'res precipitation = .TRUE.'
3398                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
3399             ENDIF
3400             unit = 'kg/kg m/s'
3401
3402          CASE ( 'q', 'vpt' )
3403             IF (  .NOT.  humidity )  THEN
3404                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3405                                 'res humidity = .TRUE.'
3406                CALL message( 'check_parameters', 'PA0105', 1, 2, 0, 6, 0 )
3407             ENDIF
3408             IF ( TRIM( var ) == 'q'   )  unit = 'kg/kg'
3409             IF ( TRIM( var ) == 'vpt' )  unit = 'K'
3410
3411          CASE ( 'qc' )
3412             IF (  .NOT.  cloud_physics )  THEN
3413                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3414                         'res cloud_physics = .TRUE.'
3415                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3416             ELSEIF ( icloud_scheme /= 0 ) THEN
3417                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3418                         'res cloud_scheme = seifert_beheng'
3419                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
3420             ENDIF
3421             unit = 'kg/kg'
3422
3423          CASE ( 'ql' )
3424             IF ( .NOT.  ( cloud_physics  .OR.  cloud_droplets ) )  THEN
3425                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3426                         'res cloud_physics = .TRUE. or cloud_droplets = .TRUE.'
3427                CALL message( 'check_parameters', 'PA0106', 1, 2, 0, 6, 0 )
3428             ENDIF
3429             unit = 'kg/kg'
3430
3431          CASE ( 'ql_c', 'ql_v', 'ql_vp' )
3432             IF (  .NOT.  cloud_droplets )  THEN
3433                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3434                                 'res cloud_droplets = .TRUE.'
3435                CALL message( 'check_parameters', 'PA0107', 1, 2, 0, 6, 0 )
3436             ENDIF
3437             IF ( TRIM( var ) == 'ql_c'  )  unit = 'kg/kg'
3438             IF ( TRIM( var ) == 'ql_v'  )  unit = 'm3'
3439             IF ( TRIM( var ) == 'ql_vp' )  unit = 'none'
3440
3441          CASE ( 'qr' )
3442             IF (  .NOT.  cloud_physics )  THEN
3443                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3444                         'res cloud_physics = .TRUE.'
3445                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3446             ELSEIF ( icloud_scheme /= 0 ) THEN
3447                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3448                         'res cloud_scheme = seifert_beheng'
3449                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
3450             ELSEIF (  .NOT.  precipitation )  THEN
3451                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3452                                 'res precipitation = .TRUE.'
3453                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
3454             ENDIF
3455             unit = 'kg/kg'
3456
3457          CASE ( 'qv' )
3458             IF (  .NOT.  cloud_physics )  THEN
3459                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3460                                 'res cloud_physics = .TRUE.'
3461                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3462             ENDIF
3463             unit = 'kg/kg'
3464
3465          CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_lw_cs_hr', 'rad_lw_hr',       &
3466                 'rad_sw_in', 'rad_sw_out', 'rad_sw_cs_hr', 'rad_sw_hr' )
3467             IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
3468                message_string = '"output of "' // TRIM( var ) // '" requi' // &
3469                                 'res radiation = .TRUE. and ' //              &
3470                                 'radiation_scheme = "rrtmg"'
3471                CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
3472             ENDIF
3473             unit = 'W/m2'
3474
3475          CASE ( 'rho' )
3476             IF (  .NOT.  ocean )  THEN
3477                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3478                                 'res ocean = .TRUE.'
3479                CALL message( 'check_parameters', 'PA0109', 1, 2, 0, 6, 0 )
3480             ENDIF
3481             unit = 'kg/m3'
3482
3483          CASE ( 's' )
3484             IF (  .NOT.  passive_scalar )  THEN
3485                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3486                                 'res passive_scalar = .TRUE.'
3487                CALL message( 'check_parameters', 'PA0110', 1, 2, 0, 6, 0 )
3488             ENDIF
3489             unit = 'conc'
3490
3491          CASE ( 'sa' )
3492             IF (  .NOT.  ocean )  THEN
3493                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3494                                 'res ocean = .TRUE.'
3495                CALL message( 'check_parameters', 'PA0109', 1, 2, 0, 6, 0 )
3496             ENDIF
3497             unit = 'psu'
3498
3499          CASE ( 't_soil' )
3500             IF (  .NOT.  land_surface )  THEN
3501                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3502                         'land_surface = .TRUE.'
3503                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3504             ENDIF
3505             unit = 'K'
3506
3507
3508          CASE ( 'c_liq*', 'c_soil*', 'c_veg*', 'ghf_eb*', 'lai*', 'lwp*',     &
3509                 'm_liq_eb*', 'ol*', 'pra*', 'prr*', 'qsws*', 'qsws_eb*',      &
3510                 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*', 'rad_net*',  &
3511                 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*', 'rrtm_asdir*',   &
3512                 'r_a*', 'r_s*', 'shf*', 'shf_eb*', 't*', 'u*', 'z0*', 'z0h*', &
3513                 'z0q*' )
3514             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
3515                message_string = 'illegal value for data_output: "' //         &
3516                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
3517                                 'cross sections are allowed for this value'
3518                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
3519             ENDIF
3520             IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
3521                IF ( TRIM( var ) == 'rrtm_aldif*'  .OR.                        &
3522                     TRIM( var ) == 'rrtm_aldir*'  .OR.                        &
3523                     TRIM( var ) == 'rrtm_asdif*'  .OR.                        &
3524                     TRIM( var ) == 'rrtm_asdir*'      )                       &
3525                THEN
3526                   message_string = 'output of "' // TRIM( var ) // '" require'&
3527                                    // 's radiation = .TRUE. and radiation_sch'&
3528                                    // 'eme = "rrtmg"'
3529                   CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
3530                ENDIF
3531             ENDIF
3532
3533             IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT.  land_surface )  THEN
3534                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3535                                 'res land_surface = .TRUE.'
3536                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3537             ENDIF
3538             IF ( TRIM( var ) == 'c_soil*'  .AND.  .NOT.  land_surface )  THEN
3539                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3540                                 'res land_surface = .TRUE.'
3541                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3542             ENDIF
3543             IF ( TRIM( var ) == 'c_veg*'  .AND.  .NOT. land_surface )  THEN
3544                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3545                                 'res land_surface = .TRUE.'
3546                CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
3547             ENDIF
3548             IF ( TRIM( var ) == 'ghf_eb*'  .AND.  .NOT.  land_surface )  THEN
3549                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3550                                 'res land_surface = .TRUE.'
3551                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3552             ENDIF
3553             IF ( TRIM( var ) == 'lai*'  .AND.  .NOT.  land_surface )  THEN
3554                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3555                                 'res land_surface = .TRUE.'
3556                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3557             ENDIF
3558             IF ( TRIM( var ) == 'lwp*'  .AND.  .NOT. cloud_physics )  THEN
3559                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3560                                 'res cloud_physics = .TRUE.'
3561                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3562             ENDIF
3563             IF ( TRIM( var ) == 'm_liq_eb*'  .AND.  .NOT.  land_surface )  THEN
3564                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3565                                 'res land_surface = .TRUE.'
3566                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3567             ENDIF
3568             IF ( TRIM( var ) == 'pra*'  .AND.  .NOT. precipitation )  THEN
3569                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3570                                 'res precipitation = .TRUE.'
3571                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
3572             ENDIF
3573             IF ( TRIM( var ) == 'pra*'  .AND.  j == 1 )  THEN
3574                message_string = 'temporal averaging of precipitation ' //     &
3575                          'amount "' // TRIM( var ) // '" is not possible'
3576                CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 )
3577             ENDIF
3578             IF ( TRIM( var ) == 'prr*'  .AND.  .NOT.  precipitation )  THEN
3579                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3580                                 'res precipitation = .TRUE.'
3581                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
3582             ENDIF
3583             IF ( TRIM( var ) == 'qsws*'  .AND.  .NOT.  humidity )  THEN
3584                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3585                                 'res humidity = .TRUE.'
3586                CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 )
3587             ENDIF
3588             IF ( TRIM( var ) == 'qsws_eb*'  .AND.  .NOT.  land_surface )  THEN
3589                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3590                                 'res land_surface = .TRUE.'
3591                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3592             ENDIF
3593             IF ( TRIM( var ) == 'qsws_liq_eb*'  .AND.  .NOT. land_surface )   &
3594             THEN
3595                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3596                                 'res land_surface = .TRUE.'
3597                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3598             ENDIF
3599             IF ( TRIM( var ) == 'qsws_soil_eb*'  .AND.  .NOT.  land_surface ) &
3600             THEN
3601                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3602                                 'res land_surface = .TRUE.'
3603                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3604             ENDIF
3605             IF ( TRIM( var ) == 'qsws_veg_eb*'  .AND.  .NOT. land_surface )   &
3606             THEN
3607                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3608                                 'res land_surface = .TRUE.'
3609                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3610             ENDIF
3611             IF ( TRIM( var ) == 'r_a*'  .AND.  .NOT.  land_surface ) &
3612             THEN
3613                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3614                                 'res land_surface = .TRUE.'
3615                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3616             ENDIF
3617             IF ( TRIM( var ) == 'r_s*'  .AND.  .NOT.  land_surface ) &
3618             THEN
3619                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3620                                 'res land_surface = .TRUE.'
3621                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3622             ENDIF
3623
3624             IF ( TRIM( var ) == 'c_liq*' )  unit = 'none'
3625             IF ( TRIM( var ) == 'c_soil*')  unit = 'none'
3626             IF ( TRIM( var ) == 'c_veg*' )  unit = 'none'
3627             IF ( TRIM( var ) == 'ghf_eb*')  unit = 'W/m2'
3628             IF ( TRIM( var ) == 'lai*'   )  unit = 'none'
3629             IF ( TRIM( var ) == 'lwp*'   )  unit = 'kg/m2'
3630             IF ( TRIM( var ) == 'm_liq_eb*'     )  unit = 'm'
3631             IF ( TRIM( var ) == 'ol*'   )   unit = 'm'
3632             IF ( TRIM( var ) == 'pra*'   )  unit = 'mm'
3633             IF ( TRIM( var ) == 'prr*'   )  unit = 'mm/s'
3634             IF ( TRIM( var ) == 'qsws*'  )  unit = 'kgm/kgs'
3635             IF ( TRIM( var ) == 'qsws_eb*'      ) unit = 'W/m2'
3636             IF ( TRIM( var ) == 'qsws_liq_eb*'  ) unit = 'W/m2'
3637             IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2'
3638             IF ( TRIM( var ) == 'qsws_veg_eb*'  ) unit = 'W/m2'
3639             IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'   
3640             IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''   
3641             IF ( TRIM( var ) == 'rrtm_aldir*'   ) unit = '' 
3642             IF ( TRIM( var ) == 'rrtm_asdif*'   ) unit = '' 
3643             IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = '' 
3644             IF ( TRIM( var ) == 'r_a*')     unit = 's/m'     
3645             IF ( TRIM( var ) == 'r_s*')     unit = 's/m' 
3646             IF ( TRIM( var ) == 'shf*'   )  unit = 'K*m/s'
3647             IF ( TRIM( var ) == 'shf_eb*')  unit = 'W/m2'
3648             IF ( TRIM( var ) == 't*'     )  unit = 'K'
3649             IF ( TRIM( var ) == 'u*'     )  unit = 'm/s'
3650             IF ( TRIM( var ) == 'z0*'    )  unit = 'm'
3651             IF ( TRIM( var ) == 'z0h*'   )  unit = 'm'
3652
3653
3654          CASE ( 'p', 'pt', 'u', 'v', 'w' )
3655             IF ( TRIM( var ) == 'p'  )  unit = 'Pa'
3656             IF ( TRIM( var ) == 'pt' )  unit = 'K'
3657             IF ( TRIM( var ) == 'u'  )  unit = 'm/s'
3658             IF ( TRIM( var ) == 'v'  )  unit = 'm/s'
3659             IF ( TRIM( var ) == 'w'  )  unit = 'm/s'
3660             CONTINUE
3661
3662          CASE DEFAULT
3663
3664             CALL user_check_data_output( var, unit )
3665
3666             IF ( unit == 'illegal' )  THEN
3667                IF ( data_output_user(1) /= ' ' )  THEN
3668                   message_string = 'illegal value for data_output or ' //     &
3669                         'data_output_user = "' // TRIM( data_output(i) ) // '"'
3670                   CALL message( 'check_parameters', 'PA0114', 1, 2, 0, 6, 0 )
3671                ELSE
3672                   message_string = 'illegal value for data_output = "' //     &
3673                                    TRIM( data_output(i) ) // '"'
3674                   CALL message( 'check_parameters', 'PA0115', 1, 2, 0, 6, 0 )
3675                ENDIF
3676             ENDIF
3677
3678       END SELECT
3679!
3680!--    Set the internal steering parameters appropriately
3681       IF ( k == 0 )  THEN
3682          do3d_no(j)              = do3d_no(j) + 1
3683          do3d(j,do3d_no(j))      = data_output(i)
3684          do3d_unit(j,do3d_no(j)) = unit
3685       ELSE
3686          do2d_no(j)              = do2d_no(j) + 1
3687          do2d(j,do2d_no(j))      = data_output(i)
3688          do2d_unit(j,do2d_no(j)) = unit
3689          IF ( data_output(i)(ilen-2:ilen) == '_xy' )  THEN
3690             data_output_xy(j) = .TRUE.
3691          ENDIF
3692          IF ( data_output(i)(ilen-2:ilen) == '_xz' )  THEN
3693             data_output_xz(j) = .TRUE.
3694          ENDIF
3695          IF ( data_output(i)(ilen-2:ilen) == '_yz' )  THEN
3696             data_output_yz(j) = .TRUE.
3697          ENDIF
3698       ENDIF
3699
3700       IF ( j == 1 )  THEN
3701!
3702!--       Check, if variable is already subject to averaging
3703          found = .FALSE.
3704          DO  k = 1, doav_n
3705             IF ( TRIM( doav(k) ) == TRIM( var ) )  found = .TRUE.
3706          ENDDO
3707
3708          IF ( .NOT. found )  THEN
3709             doav_n = doav_n + 1
3710             doav(doav_n) = var
3711          ENDIF
3712       ENDIF
3713
3714       i = i + 1
3715    ENDDO
3716
3717!
3718!-- Averaged 2d or 3d output requires that an averaging interval has been set
3719    IF ( doav_n > 0  .AND.  averaging_interval == 0.0_wp )  THEN
3720       WRITE( message_string, * )  'output of averaged quantity "',            &
3721                                   TRIM( doav(1) ), '_av" requires to set a ', &
3722                                   'non-zero & averaging interval'
3723       CALL message( 'check_parameters', 'PA0323', 1, 2, 0, 6, 0 )
3724    ENDIF
3725
3726!
3727!-- Check sectional planes and store them in one shared array
3728    IF ( ANY( section_xy > nz + 1 ) )  THEN
3729       WRITE( message_string, * )  'section_xy must be <= nz + 1 = ', nz + 1
3730       CALL message( 'check_parameters', 'PA0319', 1, 2, 0, 6, 0 )
3731    ENDIF
3732    IF ( ANY( section_xz > ny + 1 ) )  THEN
3733       WRITE( message_string, * )  'section_xz must be <= ny + 1 = ', ny + 1
3734       CALL message( 'check_parameters', 'PA0320', 1, 2, 0, 6, 0 )
3735    ENDIF
3736    IF ( ANY( section_yz > nx + 1 ) )  THEN
3737       WRITE( message_string, * )  'section_yz must be <= nx + 1 = ', nx + 1
3738       CALL message( 'check_parameters', 'PA0321', 1, 2, 0, 6, 0 )
3739    ENDIF
3740    section(:,1) = section_xy
3741    section(:,2) = section_xz
3742    section(:,3) = section_yz
3743
3744!
3745!-- Upper plot limit for 2D vertical sections
3746    IF ( z_max_do2d == -1.0_wp )  z_max_do2d = zu(nzt)
3747    IF ( z_max_do2d < zu(nzb+1)  .OR.  z_max_do2d > zu(nzt) )  THEN
3748       WRITE( message_string, * )  'z_max_do2d = ', z_max_do2d,                &
3749                    ' must be >= ', zu(nzb+1), '(zu(nzb+1)) and <= ', zu(nzt), &
3750                    ' (zu(nzt))'
3751       CALL message( 'check_parameters', 'PA0116', 1, 2, 0, 6, 0 )
3752    ENDIF
3753
3754!
3755!-- Upper plot limit for 3D arrays
3756    IF ( nz_do3d == -9999 )  nz_do3d = nzt + 1
3757
3758!
3759!-- Set output format string (used in header)
3760    SELECT CASE ( netcdf_data_format )
3761       CASE ( 1 )
3762          netcdf_data_format_string = 'netCDF classic'
3763       CASE ( 2 )
3764          netcdf_data_format_string = 'netCDF 64bit offset'
3765       CASE ( 3 )
3766          netcdf_data_format_string = 'netCDF4/HDF5'
3767       CASE ( 4 )
3768          netcdf_data_format_string = 'netCDF4/HDF5 classic'
3769       CASE ( 5 )
3770          netcdf_data_format_string = 'parallel netCDF4/HDF5'
3771       CASE ( 6 )
3772          netcdf_data_format_string = 'parallel netCDF4/HDF5 classic'
3773
3774    END SELECT
3775
3776!
3777!-- Check mask conditions
3778    DO mid = 1, max_masks
3779       IF ( data_output_masks(mid,1) /= ' '  .OR.                              &
3780            data_output_masks_user(mid,1) /= ' ' ) THEN
3781          masks = masks + 1
3782       ENDIF
3783    ENDDO
3784   
3785    IF ( masks < 0  .OR.  masks > max_masks )  THEN
3786       WRITE( message_string, * )  'illegal value: masks must be >= 0 and ',   &
3787            '<= ', max_masks, ' (=max_masks)'
3788       CALL message( 'check_parameters', 'PA0325', 1, 2, 0, 6, 0 )
3789    ENDIF
3790    IF ( masks > 0 )  THEN
3791       mask_scale(1) = mask_scale_x
3792       mask_scale(2) = mask_scale_y
3793       mask_scale(3) = mask_scale_z
3794       IF ( ANY( mask_scale <= 0.0_wp ) )  THEN
3795          WRITE( message_string, * )                                           &
3796               'illegal value: mask_scale_x, mask_scale_y and mask_scale_z',   &
3797               'must be > 0.0'
3798          CALL message( 'check_parameters', 'PA0326', 1, 2, 0, 6, 0 )
3799       ENDIF
3800!
3801!--    Generate masks for masked data output
3802!--    Parallel netcdf output is not tested so far for masked data, hence
3803!--    netcdf_data_format is switched back to non-paralell output.
3804       netcdf_data_format_save = netcdf_data_format
3805       IF ( netcdf_data_format > 4 )  THEN
3806          IF ( netcdf_data_format == 5 ) netcdf_data_format = 3
3807          IF ( netcdf_data_format == 6 ) netcdf_data_format = 4
3808          message_string = 'netCDF file formats '//                            &
3809                           '5 (parallel netCDF 4) and ' //                     &
3810                           '6 (parallel netCDF 4 Classic model) '//            &
3811                           '&are currently not supported (not yet tested) ' // &
3812                           'for masked data.&Using respective non-parallel' // & 
3813                           ' output for masked data.'
3814          CALL message( 'check_parameters', 'PA0383', 0, 0, 0, 6, 0 )
3815       ENDIF
3816       CALL init_masks
3817       netcdf_data_format = netcdf_data_format_save
3818    ENDIF
3819
3820!
3821!-- Check the NetCDF data format
3822#if ! defined ( __check )
3823    IF ( netcdf_data_format > 2 )  THEN
3824#if defined( __netcdf4 )
3825       CONTINUE
3826#else
3827       message_string = 'netCDF: netCDF4 format requested but no ' //          &
3828                        'cpp-directive __netcdf4 given & switch '  //          &
3829                        'back to 64-bit offset format'
3830       CALL message( 'check_parameters', 'PA0171', 0, 1, 0, 6, 0 )
3831       netcdf_data_format = 2
3832#endif
3833    ENDIF
3834    IF ( netcdf_data_format > 4 )  THEN
3835#if defined( __netcdf4 ) && defined( __netcdf4_parallel )
3836       CONTINUE
3837#else
3838       message_string = 'netCDF: netCDF4 parallel output requested but no ' // &
3839                        'cpp-directive __netcdf4_parallel given & switch '  // &
3840                        'back to netCDF4 non-parallel output'
3841       CALL message( 'check_parameters', 'PA0099', 0, 1, 0, 6, 0 )
3842       netcdf_data_format = netcdf_data_format - 2
3843#endif
3844    ENDIF
3845#endif
3846
3847!
3848!-- Calculate fixed number of output time levels for parallel netcdf output.
3849!-- The time dimension has to be defined as limited for parallel output,
3850!-- because otherwise the I/O performance drops significantly.
3851    IF ( netcdf_data_format > 4 )  THEN
3852
3853!
3854!--    Check if any of the follwoing data output interval is 0.0s, which is
3855!--    not allowed for parallel output.
3856       CALL check_dt_do( dt_do3d,    'dt_do3d'    )
3857       CALL check_dt_do( dt_do2d_xy, 'dt_do2d_xy' )
3858       CALL check_dt_do( dt_do2d_xz, 'dt_do2d_xz' )
3859       CALL check_dt_do( dt_do2d_yz, 'dt_do2d_yz' )
3860
3861       ntdim_3d(0)    = INT( ( end_time - skip_time_do3d ) / dt_do3d )
3862       IF ( do3d_at_begin ) ntdim_3d(0) = ntdim_3d(0) + 1
3863       ntdim_3d(1)    = INT( ( end_time - skip_time_data_output_av )           &
3864                             / dt_data_output_av )
3865       ntdim_2d_xy(0) = INT( ( end_time - skip_time_do2d_xy ) / dt_do2d_xy )
3866       ntdim_2d_xz(0) = INT( ( end_time - skip_time_do2d_xz ) / dt_do2d_xz )
3867       ntdim_2d_yz(0) = INT( ( end_time - skip_time_do2d_yz ) / dt_do2d_yz )
3868       IF ( do2d_at_begin )  THEN
3869          ntdim_2d_xy(0) = ntdim_2d_xy(0) + 1
3870          ntdim_2d_xz(0) = ntdim_2d_xz(0) + 1
3871          ntdim_2d_yz(0) = ntdim_2d_yz(0) + 1
3872       ENDIF
3873       ntdim_2d_xy(1) = ntdim_3d(1)
3874       ntdim_2d_xz(1) = ntdim_3d(1)
3875       ntdim_2d_yz(1) = ntdim_3d(1)
3876
3877    ENDIF
3878
3879!
3880!-- Check, whether a constant diffusion coefficient shall be used
3881    IF ( km_constant /= -1.0_wp )  THEN
3882       IF ( km_constant < 0.0_wp )  THEN
3883          WRITE( message_string, * )  'km_constant = ', km_constant, ' < 0.0'
3884          CALL message( 'check_parameters', 'PA0121', 1, 2, 0, 6, 0 )
3885       ELSE
3886          IF ( prandtl_number < 0.0_wp )  THEN
3887             WRITE( message_string, * )  'prandtl_number = ', prandtl_number,  &
3888                                         ' < 0.0'
3889             CALL message( 'check_parameters', 'PA0122', 1, 2, 0, 6, 0 )
3890          ENDIF
3891          constant_diffusion = .TRUE.
3892
3893          IF ( constant_flux_layer )  THEN
3894             message_string = 'constant_flux_layer is not allowed with fixed ' &
3895                              // 'value of km'
3896             CALL message( 'check_parameters', 'PA0123', 1, 2, 0, 6, 0 )
3897          ENDIF
3898       ENDIF
3899    ENDIF
3900
3901!
3902!-- In case of non-cyclic lateral boundaries and a damping layer for the
3903!-- potential temperature, check the width of the damping layer
3904    IF ( bc_lr /= 'cyclic' ) THEN
3905       IF ( pt_damping_width < 0.0_wp  .OR.                                    &
3906            pt_damping_width > REAL( nx * dx ) )  THEN
3907          message_string = 'pt_damping_width out of range'
3908          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
3909       ENDIF
3910    ENDIF
3911
3912    IF ( bc_ns /= 'cyclic' )  THEN
3913       IF ( pt_damping_width < 0.0_wp  .OR.                                    &
3914            pt_damping_width > REAL( ny * dy ) )  THEN
3915          message_string = 'pt_damping_width out of range'
3916          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
3917       ENDIF
3918    ENDIF
3919
3920!
3921!-- Check value range for zeta = z/L
3922    IF ( zeta_min >= zeta_max )  THEN
3923       WRITE( message_string, * )  'zeta_min = ', zeta_min, ' must be less ',  &
3924                                   'than zeta_max = ', zeta_max
3925       CALL message( 'check_parameters', 'PA0125', 1, 2, 0, 6, 0 )
3926    ENDIF
3927
3928!
3929!-- Check random generator
3930    IF ( (random_generator /= 'system-specific'      .AND.                     &
3931          random_generator /= 'random-parallel'   )  .AND.                     &
3932          random_generator /= 'numerical-recipes' )  THEN
3933       message_string = 'unknown random generator: random_generator = "' //    &
3934                        TRIM( random_generator ) // '"'
3935       CALL message( 'check_parameters', 'PA0135', 1, 2, 0, 6, 0 )
3936    ENDIF
3937
3938!
3939!-- Determine upper and lower hight level indices for random perturbations
3940    IF ( disturbance_level_b == -9999999.9_wp )  THEN
3941       IF ( ocean )  THEN
3942          disturbance_level_b     = zu((nzt*2)/3)
3943          disturbance_level_ind_b = ( nzt * 2 ) / 3
3944       ELSE
3945          disturbance_level_b     = zu(nzb+3)
3946          disturbance_level_ind_b = nzb + 3
3947       ENDIF
3948    ELSEIF ( disturbance_level_b < zu(3) )  THEN
3949       WRITE( message_string, * )  'disturbance_level_b = ',                   &
3950                           disturbance_level_b, ' must be >= ', zu(3), '(zu(3))'
3951       CALL message( 'check_parameters', 'PA0126', 1, 2, 0, 6, 0 )
3952    ELSEIF ( disturbance_level_b > zu(nzt-2) )  THEN
3953       WRITE( message_string, * )  'disturbance_level_b = ',                   &
3954                   disturbance_level_b, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
3955       CALL message( 'check_parameters', 'PA0127', 1, 2, 0, 6, 0 )
3956    ELSE
3957       DO  k = 3, nzt-2
3958          IF ( disturbance_level_b <= zu(k) )  THEN
3959             disturbance_level_ind_b = k
3960             EXIT
3961          ENDIF
3962       ENDDO
3963    ENDIF
3964
3965    IF ( disturbance_level_t == -9999999.9_wp )  THEN
3966       IF ( ocean )  THEN
3967          disturbance_level_t     = zu(nzt-3)
3968          disturbance_level_ind_t = nzt - 3
3969       ELSE
3970          disturbance_level_t     = zu(nzt/3)
3971          disturbance_level_ind_t = nzt / 3
3972       ENDIF
3973    ELSEIF ( disturbance_level_t > zu(nzt-2) )  THEN
3974       WRITE( message_string, * )  'disturbance_level_t = ',                   &
3975                   disturbance_level_t, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
3976       CALL message( 'check_parameters', 'PA0128', 1, 2, 0, 6, 0 )
3977    ELSEIF ( disturbance_level_t < disturbance_level_b )  THEN
3978       WRITE( message_string, * )  'disturbance_level_t = ',                   &
3979                   disturbance_level_t, ' must be >= disturbance_level_b = ',  &
3980                   disturbance_level_b
3981       CALL message( 'check_parameters', 'PA0129', 1, 2, 0, 6, 0 )
3982    ELSE
3983       DO  k = 3, nzt-2
3984          IF ( disturbance_level_t <= zu(k) )  THEN
3985             disturbance_level_ind_t = k
3986             EXIT
3987          ENDIF
3988       ENDDO
3989    ENDIF
3990
3991!
3992!-- Check again whether the levels determined this way are ok.
3993!-- Error may occur at automatic determination and too few grid points in
3994!-- z-direction.
3995    IF ( disturbance_level_ind_t < disturbance_level_ind_b )  THEN
3996       WRITE( message_string, * )  'disturbance_level_ind_t = ',               &
3997                disturbance_level_ind_t, ' must be >= disturbance_level_b = ', &
3998                disturbance_level_b
3999       CALL message( 'check_parameters', 'PA0130', 1, 2, 0, 6, 0 )
4000    ENDIF
4001
4002!
4003!-- Determine the horizontal index range for random perturbations.
4004!-- In case of non-cyclic horizontal boundaries, no perturbations are imposed
4005!-- near the inflow and the perturbation area is further limited to ...(1)
4006!-- after the initial phase of the flow.
4007   
4008    IF ( bc_lr /= 'cyclic' )  THEN
4009       IF ( inflow_disturbance_begin == -1 )  THEN
4010          inflow_disturbance_begin = MIN( 10, nx/2 )
4011       ENDIF
4012       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > nx )&
4013       THEN
4014          message_string = 'inflow_disturbance_begin out of range'
4015          CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
4016       ENDIF
4017       IF ( inflow_disturbance_end == -1 )  THEN
4018          inflow_disturbance_end = MIN( 100, 3*nx/4 )
4019       ENDIF
4020       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > nx )    &
4021       THEN
4022          message_string = 'inflow_disturbance_end out of range'
4023          CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
4024       ENDIF
4025    ELSEIF ( bc_ns /= 'cyclic' )  THEN
4026       IF ( inflow_disturbance_begin == -1 )  THEN
4027          inflow_disturbance_begin = MIN( 10, ny/2 )
4028       ENDIF
4029       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > ny )&
4030       THEN
4031          message_string = 'inflow_disturbance_begin out of range'
4032          CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
4033       ENDIF
4034       IF ( inflow_disturbance_end == -1 )  THEN
4035          inflow_disturbance_end = MIN( 100, 3*ny/4 )
4036       ENDIF
4037       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > ny )    &
4038       THEN
4039          message_string = 'inflow_disturbance_end out of range'
4040          CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
4041       ENDIF
4042    ENDIF
4043
4044    IF ( random_generator == 'random-parallel' )  THEN
4045       dist_nxl = nxl;  dist_nxr = nxr
4046       dist_nys = nys;  dist_nyn = nyn
4047       IF ( bc_lr == 'radiation/dirichlet' )  THEN
4048          dist_nxr    = MIN( nx - inflow_disturbance_begin, nxr )
4049          dist_nxl(1) = MAX( nx - inflow_disturbance_end, nxl )
4050       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
4051          dist_nxl    = MAX( inflow_disturbance_begin, nxl )
4052          dist_nxr(1) = MIN( inflow_disturbance_end, nxr )
4053       ENDIF
4054       IF ( bc_ns == 'dirichlet/radiation' )  THEN
4055          dist_nyn    = MIN( ny - inflow_disturbance_begin, nyn )
4056          dist_nys(1) = MAX( ny - inflow_disturbance_end, nys )
4057       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
4058          dist_nys    = MAX( inflow_disturbance_begin, nys )
4059          dist_nyn(1) = MIN( inflow_disturbance_end, nyn )
4060       ENDIF
4061    ELSE
4062       dist_nxl = 0;  dist_nxr = nx
4063       dist_nys = 0;  dist_nyn = ny
4064       IF ( bc_lr == 'radiation/dirichlet' )  THEN
4065          dist_nxr    = nx - inflow_disturbance_begin
4066          dist_nxl(1) = nx - inflow_disturbance_end
4067       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
4068          dist_nxl    = inflow_disturbance_begin
4069          dist_nxr(1) = inflow_disturbance_end
4070       ENDIF
4071       IF ( bc_ns == 'dirichlet/radiation' )  THEN
4072          dist_nyn    = ny - inflow_disturbance_begin
4073          dist_nys(1) = ny - inflow_disturbance_end
4074       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
4075          dist_nys    = inflow_disturbance_begin
4076          dist_nyn(1) = inflow_disturbance_end
4077       ENDIF
4078    ENDIF
4079
4080!
4081!-- A turbulent inflow requires Dirichlet conditions at the respective inflow
4082!-- boundary (so far, a turbulent inflow is realized from the left side only)
4083    IF ( turbulent_inflow  .AND.  bc_lr /= 'dirichlet/radiation' )  THEN
4084       message_string = 'turbulent_inflow = .T. requires a Dirichlet ' //      &
4085                        'condition at the inflow boundary'
4086       CALL message( 'check_parameters', 'PA0133', 1, 2, 0, 6, 0 )
4087    ENDIF
4088
4089!
4090!-- Turbulent inflow requires that 3d arrays have been cyclically filled with
4091!-- data from prerun in the first main run
4092    IF ( turbulent_inflow  .AND.  initializing_actions /= 'cyclic_fill'  .AND. &
4093         initializing_actions /= 'read_restart_data' )  THEN
4094       message_string = 'turbulent_inflow = .T. requires ' //                  &
4095                        'initializing_actions = ''cyclic_fill'' '
4096       CALL message( 'check_parameters', 'PA0055', 1, 2, 0, 6, 0 )
4097    ENDIF
4098
4099!
4100!-- In case of turbulent inflow calculate the index of the recycling plane
4101    IF ( turbulent_inflow )  THEN
4102       IF ( recycling_width == 9999999.9_wp )  THEN
4103!
4104!--       Set the default value for the width of the recycling domain
4105          recycling_width = 0.1_wp * nx * dx
4106       ELSE
4107          IF ( recycling_width < dx  .OR.  recycling_width > nx * dx )  THEN
4108             WRITE( message_string, * )  'illegal value for recycling_width:', &
4109                                         ' ', recycling_width
4110             CALL message( 'check_parameters', 'PA0134', 1, 2, 0, 6, 0 )
4111          ENDIF
4112       ENDIF
4113!
4114!--    Calculate the index
4115       recycling_plane = recycling_width / dx
4116    ENDIF
4117
4118!
4119!-- Determine damping level index for 1D model
4120    IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
4121       IF ( damp_level_1d == -1.0_wp )  THEN
4122          damp_level_1d     = zu(nzt+1)
4123          damp_level_ind_1d = nzt + 1
4124       ELSEIF ( damp_level_1d < 0.0_wp  .OR.  damp_level_1d > zu(nzt+1) )  THEN
4125          WRITE( message_string, * )  'damp_level_1d = ', damp_level_1d,       &
4126                 ' must be > 0.0 and < ', zu(nzt+1), '(zu(nzt+1))'
4127          CALL message( 'check_parameters', 'PA0136', 1, 2, 0, 6, 0 )
4128       ELSE
4129          DO  k = 1, nzt+1
4130             IF ( damp_level_1d <= zu(k) )  THEN
4131                damp_level_ind_1d = k
4132                EXIT
4133             ENDIF
4134          ENDDO
4135       ENDIF
4136    ENDIF
4137
4138!
4139!-- Check some other 1d-model parameters
4140    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND.                   &
4141         TRIM( mixing_length_1d ) /= 'blackadar' )  THEN
4142       message_string = 'mixing_length_1d = "' // TRIM( mixing_length_1d ) //  &
4143                        '" is unknown'
4144       CALL message( 'check_parameters', 'PA0137', 1, 2, 0, 6, 0 )
4145    ENDIF
4146    IF ( TRIM( dissipation_1d ) /= 'as_in_3d_model'  .AND.                     &
4147         TRIM( dissipation_1d ) /= 'detering' )  THEN
4148       message_string = 'dissipation_1d = "' // TRIM( dissipation_1d ) //      &
4149                        '" is unknown'
4150       CALL message( 'check_parameters', 'PA0138', 1, 2, 0, 6, 0 )
4151    ENDIF
4152
4153!
4154!-- Set time for the next user defined restart (time_restart is the
4155!-- internal parameter for steering restart events)
4156    IF ( restart_time /= 9999999.9_wp )  THEN
4157       IF ( restart_time > time_since_reference_point )  THEN
4158          time_restart = restart_time
4159       ENDIF
4160    ELSE
4161!
4162!--    In case of a restart run, set internal parameter to default (no restart)
4163!--    if the NAMELIST-parameter restart_time is omitted
4164       time_restart = 9999999.9_wp
4165    ENDIF
4166
4167!
4168!-- Set default value of the time needed to terminate a model run
4169    IF ( termination_time_needed == -1.0_wp )  THEN
4170       IF ( host(1:3) == 'ibm' )  THEN
4171          termination_time_needed = 300.0_wp
4172       ELSE
4173          termination_time_needed = 35.0_wp
4174       ENDIF
4175    ENDIF
4176
4177!
4178!-- Check the time needed to terminate a model run
4179    IF ( host(1:3) == 't3e' )  THEN
4180!
4181!--    Time needed must be at least 30 seconds on all CRAY machines, because
4182!--    MPP_TREMAIN gives the remaining CPU time only in steps of 30 seconds
4183       IF ( termination_time_needed <= 30.0_wp )  THEN
4184          WRITE( message_string, * )  'termination_time_needed = ',            &
4185                 termination_time_needed, ' must be > 30.0 on host "',         &
4186                 TRIM( host ), '"'
4187          CALL message( 'check_parameters', 'PA0139', 1, 2, 0, 6, 0 )
4188       ENDIF
4189    ELSEIF ( host(1:3) == 'ibm' )  THEN
4190!
4191!--    On IBM-regatta machines the time should be at least 300 seconds,
4192!--    because the job time consumed before executing palm (for compiling,
4193!--    copying of files, etc.) has to be regarded
4194       IF ( termination_time_needed < 300.0_wp )  THEN
4195          WRITE( message_string, * )  'termination_time_needed = ',            &
4196                 termination_time_needed, ' should be >= 300.0 on host "',     &
4197                 TRIM( host ), '"'
4198          CALL message( 'check_parameters', 'PA0140', 1, 2, 0, 6, 0 )
4199       ENDIF
4200    ENDIF
4201
4202!
4203!-- Check pressure gradient conditions
4204    IF ( dp_external  .AND.  conserve_volume_flow )  THEN
4205       WRITE( message_string, * )  'Both dp_external and conserve_volume_flo', &
4206            'w are .TRUE. but one of them must be .FALSE.'
4207       CALL message( 'check_parameters', 'PA0150', 1, 2, 0, 6, 0 )
4208    ENDIF
4209    IF ( dp_external )  THEN
4210       IF ( dp_level_b < zu(nzb)  .OR.  dp_level_b > zu(nzt) )  THEN
4211          WRITE( message_string, * )  'dp_level_b = ', dp_level_b, ' is out ', &
4212               ' of range'
4213          CALL message( 'check_parameters', 'PA0151', 1, 2, 0, 6, 0 )
4214       ENDIF
4215       IF ( .NOT. ANY( dpdxy /= 0.0_wp ) )  THEN
4216          WRITE( message_string, * )  'dp_external is .TRUE. but dpdxy is ze', &
4217               'ro, i.e. the external pressure gradient & will not be applied'
4218          CALL message( 'check_parameters', 'PA0152', 0, 1, 0, 6, 0 )
4219       ENDIF
4220    ENDIF
4221    IF ( ANY( dpdxy /= 0.0_wp )  .AND.  .NOT.  dp_external )  THEN
4222       WRITE( message_string, * )  'dpdxy is nonzero but dp_external is ',     &
4223            '.FALSE., i.e. the external pressure gradient & will not be applied'
4224       CALL message( 'check_parameters', 'PA0153', 0, 1, 0, 6, 0 )
4225    ENDIF
4226    IF ( conserve_volume_flow )  THEN
4227       IF ( TRIM( conserve_volume_flow_mode ) == 'default' )  THEN
4228
4229          conserve_volume_flow_mode = 'initial_profiles'
4230
4231       ELSEIF ( TRIM( conserve_volume_flow_mode ) /= 'initial_profiles' .AND.  &
4232            TRIM( conserve_volume_flow_mode ) /= 'inflow_profile' .AND.        &
4233            TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' )  THEN
4234          WRITE( message_string, * )  'unknown conserve_volume_flow_mode: ',   &
4235               conserve_volume_flow_mode
4236          CALL message( 'check_parameters', 'PA0154', 1, 2, 0, 6, 0 )
4237       ENDIF
4238       IF ( (bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic')  .AND.                &
4239          TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
4240          WRITE( message_string, * )  'non-cyclic boundary conditions ',       &
4241               'require  conserve_volume_flow_mode = ''initial_profiles'''
4242          CALL message( 'check_parameters', 'PA0155', 1, 2, 0, 6, 0 )
4243       ENDIF
4244       IF ( bc_lr == 'cyclic'  .AND.  bc_ns == 'cyclic'  .AND.                 &
4245            TRIM( conserve_volume_flow_mode ) == 'inflow_profile' )  THEN
4246          WRITE( message_string, * )  'cyclic boundary conditions ',           &
4247               'require conserve_volume_flow_mode = ''initial_profiles''',     &
4248               ' or ''bulk_velocity'''
4249          CALL message( 'check_parameters', 'PA0156', 1, 2, 0, 6, 0 )
4250       ENDIF
4251    ENDIF
4252    IF ( ( u_bulk /= 0.0_wp  .OR.  v_bulk /= 0.0_wp )  .AND.                   &
4253         ( .NOT. conserve_volume_flow  .OR.                                    &
4254         TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' ) )  THEN
4255       WRITE( message_string, * )  'nonzero bulk velocity requires ',          &
4256            'conserve_volume_flow = .T. and ',                                 &
4257            'conserve_volume_flow_mode = ''bulk_velocity'''
4258       CALL message( 'check_parameters', 'PA0157', 1, 2, 0, 6, 0 )
4259    ENDIF
4260
4261!
4262!-- Check particle attributes
4263    IF ( particle_color /= 'none' )  THEN
4264       IF ( particle_color /= 'absuv'  .AND.  particle_color /= 'pt*'  .AND.   &
4265            particle_color /= 'z' )  THEN
4266          message_string = 'illegal value for parameter particle_color: ' //   &
4267                           TRIM( particle_color)
4268          CALL message( 'check_parameters', 'PA0313', 1, 2, 0, 6, 0 )
4269       ELSE
4270          IF ( color_interval(2) <= color_interval(1) )  THEN
4271             message_string = 'color_interval(2) <= color_interval(1)'
4272             CALL message( 'check_parameters', 'PA0315', 1, 2, 0, 6, 0 )
4273          ENDIF
4274       ENDIF
4275    ENDIF
4276
4277    IF ( particle_dvrpsize /= 'none' )  THEN
4278       IF ( particle_dvrpsize /= 'absw' )  THEN
4279          message_string = 'illegal value for parameter particle_dvrpsize:' // &
4280                           ' ' // TRIM( particle_color)
4281          CALL message( 'check_parameters', 'PA0314', 1, 2, 0, 6, 0 )
4282       ELSE
4283          IF ( dvrpsize_interval(2) <= dvrpsize_interval(1) )  THEN
4284             message_string = 'dvrpsize_interval(2) <= dvrpsize_interval(1)'
4285             CALL message( 'check_parameters', 'PA0316', 1, 2, 0, 6, 0 )
4286          ENDIF
4287       ENDIF
4288    ENDIF
4289
4290!
4291!-- Check nudging and large scale forcing from external file
4292    IF ( nudging  .AND.  (  .NOT.  large_scale_forcing ) )  THEN
4293       message_string = 'Nudging requires large_scale_forcing = .T.. &'//      &
4294                        'Surface fluxes and geostrophic wind should be &'//    &
4295                        'prescribed in file LSF_DATA'
4296       CALL message( 'check_parameters', 'PA0374', 1, 2, 0, 6, 0 )
4297    ENDIF
4298
4299    IF ( large_scale_forcing  .AND.  ( bc_lr /= 'cyclic'  .OR.                 &
4300                                    bc_ns /= 'cyclic' ) )  THEN
4301       message_string = 'Non-cyclic lateral boundaries do not allow for &' //  &
4302                        'the usage of large scale forcing from external file.'
4303       CALL message( 'check_parameters', 'PA0375', 1, 2, 0, 6, 0 )
4304     ENDIF
4305
4306    IF ( large_scale_forcing  .AND.  (  .NOT.  humidity ) )  THEN
4307       message_string = 'The usage of large scale forcing from external &'//   & 
4308                        'file LSF_DATA requires humidity = .T..'
4309       CALL message( 'check_parameters', 'PA0376', 1, 2, 0, 6, 0 )
4310     ENDIF
4311
4312    IF ( large_scale_forcing  .AND.  topography /= 'flat' )  THEN
4313       message_string = 'The usage of large scale forcing from external &'//   & 
4314                        'file LSF_DATA is not implemented for non-flat topography'
4315       CALL message( 'check_parameters', 'PA0377', 1, 2, 0, 6, 0 )
4316    ENDIF
4317
4318    IF ( large_scale_forcing  .AND.  ocean  )  THEN
4319       message_string = 'The usage of large scale forcing from external &'//   & 
4320                        'file LSF_DATA is not implemented for ocean runs'
4321       CALL message( 'check_parameters', 'PA0378', 1, 2, 0, 6, 0 )
4322    ENDIF
4323
4324    CALL location_message( 'finished', .TRUE. )
4325
4326!
4327!-- Prevent empty time records in volume, cross-section and masked data in case of
4328!-- non-parallel netcdf-output in restart runs
4329    IF ( netcdf_data_format < 5 )  THEN
4330       IF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
4331          do3d_time_count    = 0
4332          do2d_xy_time_count = 0
4333          do2d_xz_time_count = 0
4334          do2d_yz_time_count = 0
4335          domask_time_count  = 0
4336       ENDIF
4337    ENDIF
4338
4339!
4340!-- Check for valid setting of most_method
4341    IF ( TRIM( most_method ) /= 'circular'  .AND.                              &
4342         TRIM( most_method ) /= 'newton'    .AND.                              &
4343         TRIM( most_method ) /= 'lookup' )  THEN
4344       message_string = 'most_method = "' // TRIM( most_method ) //            &
4345                        '" is unknown'
4346       CALL message( 'check_parameters', 'PA0416', 1, 2, 0, 6, 0 )
4347    ENDIF
4348
4349!
4350!-- Check &userpar parameters
4351    CALL user_check_parameters
4352
4353 CONTAINS
4354
4355!------------------------------------------------------------------------------!
4356! Description:
4357! ------------
4358!> Check the length of data output intervals. In case of parallel NetCDF output
4359!> the time levels of the output files need to be fixed. Therefore setting the
4360!> output interval to 0.0s (usually used to output each timestep) is not
4361!> possible as long as a non-fixed timestep is used.
4362!------------------------------------------------------------------------------!
4363
4364    SUBROUTINE check_dt_do( dt_do, dt_do_name )
4365
4366       IMPLICIT NONE
4367
4368       CHARACTER (LEN=*), INTENT (IN) :: dt_do_name !< parin variable name
4369
4370       REAL(wp), INTENT (INOUT)       :: dt_do      !< data output interval
4371
4372       IF ( dt_do == 0.0_wp )  THEN
4373          IF ( dt_fixed )  THEN
4374             WRITE( message_string, '(A,F9.4,A)' )  'Output at every '  //     &
4375                    'timestep is desired (' // dt_do_name // ' = 0.0).&'//     &
4376                    'Setting the output interval to the fixed timestep '//     &
4377                    'dt = ', dt, 's.'
4378             CALL message( 'check_parameters', 'PA0060', 0, 0, 0, 6, 0 )
4379             dt_do = dt
4380          ELSE
4381             message_string = dt_do_name // ' = 0.0 while using a ' //         &
4382                              'variable timestep and parallel netCDF4 ' //     &
4383                              'is not allowed.'
4384             CALL message( 'check_parameters', 'PA0081', 1, 2, 0, 6, 0 )
4385          ENDIF
4386       ENDIF
4387
4388    END SUBROUTINE check_dt_do
4389
4390 END SUBROUTINE check_parameters
Note: See TracBrowser for help on using the repository browser.